{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Features of BIDMat and Scala"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BIDMat is a multi-platform matrix library similar to R, Matlab, Julia or Numpy/Scipy. It takes full advantage of the very powerful Scala Language. Its intended primarily for machine learning, but is has a broad set of operations and datatypes and should be suitable for many other applications. BIDMat has several unique features:\n",
    "\n",
    "* Built from the ground up with GPU + CPU backends. BIDMat code is implementation independent. \n",
    "* GPU memory management uses caching, designed to support iterative algorithms.\n",
    "* Natural and extensible syntax (thanks to scala). Math operators include +,-,*,/,⊗,∙,∘\n",
    "* Probably the most complete support for matrix types: dense matrices of float32, double, int and long. Sparse matrices with single or double elements. All are available on CPU or GPU. \n",
    "* Highest performance sparse matrix operations on power-law data. \n",
    "\n",
    "BIDMat has several other state-of-the-art features:\n",
    "* Interactivity. Thanks to the Scala language, BIDMat is interactive and scriptable. \n",
    "* Massive code base thanks to Java. \n",
    "* Easy-to-use Parallelism, thanks to Scala's actor framework and parallel collection classes.  \n",
    "* Runs on JVM, extremely portable. Runs on Mac, Linux, Windows, Android.\n",
    "* Cluster-ready, leverages Hadoop, Yarn, Spark etc. \n",
    "\n",
    "BIDMat is a library that is loaded by a startup script, and a set of imports that include the default classes and functions. We include them explicitly in this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 CUDA device found, CUDA version 7.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.9833951,11878821888,12079398912)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import BIDMat.{CMat,CSMat,DMat,Dict,IDict,FMat,FND,GMat,GDMat,GIMat,GLMat,GSMat,GSDMat,\n",
    "               HMat,IMat,Image,LMat,Mat,ND,SMat,SBMat,SDMat}\n",
    "import BIDMat.MatFunctions._\n",
    "import BIDMat.SciFunctions._\n",
    "import BIDMat.Solvers._\n",
    "import BIDMat.JPlotting._\n",
    "\n",
    "Mat.checkMKL\n",
    "Mat.checkCUDA\n",
    "Mat.setInline\n",
    "if (Mat.hasCUDA > 0) GPUmem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These calls check that CPU and GPU native libs loaded correctly, and what GPUs are accessible.\n",
    "If you have a GPU and CUDA installed, GPUmem will printout the fraction of free memory, the absolute free memory and the total memory for the default GPU. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CPU and GPU matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BIDMat's matrix types are given in the table below. All are children of the \"Mat\" parent class, which allows code to be written generically. Many of BIDMach's learning algorithms will run with either single or double precision, dense or sparse input data. \n",
    "<table style=\"width:4in\" align=\"left\">\n",
    "<tr><td/><td colspan=\"2\"><b>CPU Matrices</b></td><td colspan=\"2\"><b>GPU Matrices</b></td></tr>\n",
    "<tr><td></td><td><b>Dense</b></td><td><b>Sparse</b></td><td><b>Dense</b></td><td><b>Sparse</b></td></tr>\n",
    "<tr><td><b>Float32</b></td><td>FMat</td><td>SMat</td><td>GMat</td><td>GSMat</td></tr>\n",
    "<tr><td><b>Float64</b></td><td>DMat</td><td>SDMat</td><td>GDMat</td><td>GSDMat</td></tr>\n",
    "<tr><td><b>Int32</b></td><td>IMat</td><td></td><td>GIMat</td><td></td></tr>\n",
    "<tr><td><b>Int64</b></td><td>LMat</td><td></td><td>GLMat</td><td></td></tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "   0.62465   0.89703   0.27582   0.13452   0.88863   0.32627   0.16042...\n",
       "   0.63891   0.44064   0.36484   0.18473   0.75125   0.18824   0.75513...\n",
       "   0.84862   0.13745   0.54177   0.69702   0.27047  0.040770  0.080366...\n",
       "   0.26990  0.088154   0.17225   0.95003   0.51944   0.89452   0.75103...\n",
       "  0.050579   0.15992   0.15957   0.49453   0.54010   0.32854   0.49945...\n",
       "   0.73763   0.38573   0.80398   0.14497   0.16867   0.89132   0.31194...\n",
       "   0.39953   0.52224   0.14720   0.11078   0.11138   0.95166  0.019181...\n",
       "   0.91040  0.076302   0.60944   0.40543   0.79685   0.85392   0.23852...\n",
       "        ..        ..        ..        ..        ..        ..        ..\n"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val n = 4096          // \"val\" designates a constant. n is statically typed (as in Int here), but its type is inferred.\n",
    "val a = rand(n,n)     // Create an nxn matrix (on the CPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BIDMat.FMat\n"
     ]
    }
   ],
   "source": [
    "%type a             // Most scientific funtions in BIDMat return single-precision results by default."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CPU matrix operations use Intel MKL acceleration for linear algebra, scientific and statistical functions. BIDMat includes \"tic\" and \"toc\" for timing, and \"flip\" and \"flop\" for floating point performance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The product took 1.32 seconds at 104 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(104.35759,1.317)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flip; val b = a * a; val gf=gflop\n",
    "print(\"The product took %4.2f seconds at %3.0f gflops\" format (gf._2, gf._1))\n",
    "gf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "GPU matrices behave very similarly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "   0.88383  0.027557   0.64643   0.67664   0.29091   0.60243   0.57059...\n",
       "   0.81384   0.76326   0.57175   0.17944   0.11371   0.22653  0.088297...\n",
       "   0.39205   0.46169   0.76375   0.85130   0.46558   0.56107   0.79527...\n",
       "   0.11010   0.51376   0.27728   0.91142  0.023597   0.16610   0.30512...\n",
       "   0.25042   0.69971   0.83499   0.46664   0.49516   0.57655   0.79105...\n",
       "   0.48823   0.72891   0.98704   0.11257   0.99232   0.10887   0.19260...\n",
       "   0.70476   0.50308   0.38125   0.80248   0.52444   0.42028   0.33877...\n",
       "   0.18731   0.11773   0.84390   0.67353   0.73904   0.19960   0.32721...\n",
       "        ..        ..        ..        ..        ..        ..        ..\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val ga = grand(n,n)            // Another nxn random matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The product took 0.04 seconds at 3124 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(3123.6125,0.044)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flip; val gb = ga * ga; val gf=gflop\n",
    "print(\"The product took %4.2f seconds at %3.0f gflops\" format (gf._2, gf._1))\n",
    "gf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BIDMat.GMat\n"
     ]
    }
   ],
   "source": [
    "%type ga"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "But much of the power of BIDMat is that we dont have to worry about matrix types. Lets explore that with an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SVD (Singular Value Decomposition) on a Budget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets try solving a real problem with this infrastructure: An approximate Singular-Value Decomposition (SVD) or PCA of a matrix $M$. We'll do this by computing the leading eigenvalues and eigenvectors of $MM^T$. The method we use is **subspace iteration** and it generalizes the *power method* for computing the largest-magnitude eigenvalue. An eigenvector is a vector $v$ such that \n",
    "$$Mv =\\lambda v$$\n",
    "where $\\lambda$ is a scalar called the eigenvalue. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "def SVD(M:Mat, ndims:Int, niter:Int) = {\n",
    "    var Q = M.zeros(M.nrows, ndims)                 // A block of ndims column vectors\n",
    "    normrnd(0, 1, Q)                                // randomly initialize the vectors\n",
    "    Mat.useCache = true                             // Turn matrix caching on\n",
    "    \n",
    "    for (i <- 0 until niter) {                      // Perform subspace iteration\n",
    "        val P = (Q.t * M *^ M).t                    // Compute P = M * M^t * Q efficiently\n",
    "        QRdecompt(P, Q, null)                       // QR-decomposition of P, saving Q\n",
    "    }\n",
    "    Mat.useCache = false                            // Turn caching off after the iteration\n",
    "    val P = (Q.t * M *^ M).t                        // Compute P again.\n",
    "    (Q, P ∙ Q)                          // Return Left singular vectors and singular values\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the code above used only the \"Mat\" matrix type. If you examine the variables V and P in a Scala IDE (Eclipse has one) you will find that they both also have type \"Mat\". Let's try it with an FMat (CPU single precision, dense matrix)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Movie Data Example\n",
    "\n",
    "We load some data from the MovieLens project. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculation took 15.66 seconds at 42.3 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  1.1839e+06\n",
       "  2.0694e+05\n",
       "  1.5514e+05\n",
       "       95136\n",
       "       91670\n",
       "       69232\n",
       "       68042\n",
       "       56050\n",
       "          ..\n"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val ndims = 32                             // Number of PCA dimension\n",
    "val niter = 128                            // Number of iterations to do\n",
    "\n",
    "val S = loadSMat(\"../data/movielens/train.smat.lz4\")(0->10000,0->4000)\n",
    "val M = full(S)                            // Put in a dense matrix\n",
    "\n",
    "flip; \n",
    "val (svecs, svals) = SVD(M, ndims, niter); // Compute the singular vectors and values\n",
    "val gf=gflop\n",
    "\n",
    "print(\"The calculation took %4.2f seconds at %2.1f gflops\" format (gf._2, gf._1))\n",
    "svals.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a peek at the singular values on a plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "450523"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S.nnz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGQCAYAAABYs5LGAAA3iUlEQVR42u2db3BV5Z3HO7Ozs7Ozsy/6pi92dnZ29t1Opzs7s7P7Zrs7GZZMGBz6R8dALFKWQIpRkKgIJaIUWzU2NikIgrZBIJumURLiH5SS0i2oUI0EKhKxwkahSEhMBSUFifw233O5t/cecv+de25yzzmfz8yZ5Nxz830uz/Ncvr/f8zt/vmAAAAAQeL5AFwAAAGDoAAAAgKEDAAAAhg4AAAAYOgAAAIYOALEv0Re+kPM20fujxDe/+c3I/tsBMHQADD3wnD171pqamiIdzABg6AAhM3T6B0MHwNABAmBYgKEDYOgAETD0dO+/cuWKrVmzxv72b//W/vzP/9z+8i//0srKypzf0+mn08r0mdzHPv/8c3vooYfs7//+7+3P/uzPEu976623bOHChfZ3f/d3ic/z5S9/2e677z775JNPMHQADB0AQ5/o/RUVFXkv2/th6F/96leve+/TTz/tmHu6zyBjHxkZYTUDAEMHCK+hezkprrW19TqTHR4etnPnzhXd0N1bX19fyv6KFSucLH79+vXXvY6hA2DoABh60mvu7HzPnj2ezNnL37S1taUsoc+dOzfl+OXLl53XZerJr//DP/wDhg6AoQNg6MmvqT7trmtPlqG7+eIXv5jTv0V1dQwdAEMHCK2he3m/VwMuhqFnqp0XasgYOgCGDhBqQ3dn6FNp6H/913+dclxn309VXwEAhg4QKEP/t3/7t5TXjh8/npO++5K2V155xXp7ewsy9Jtuuinl+LZt2657z8WLF51L6jB0AAwdAENPer+uP09+Tdd/q47+0ksvZdT/x3/8R98udYvz+uuvpxz/0pe+ZO3t7U6mrs+0d+9e57K1XA2Zu+gBYOgAkTF0XdOtG7hkM73kG7+ItWvXXvcenX1eiKELXUb3F3/xF76YMIYOgKEDRMbQxenTp23evHn2V3/1V85S+n/+53/avn37Ut6rWruburo653UZ8IIFC5zgoFBDFwMDA7Z8+XKnHCB9BRNq4ytf+YrV1NQ4mTyGDoChA0AOvP322ymm56VuDQAYOgBMIrqH+4YNG5y7w4n+/n4nS082dC2FAwBg6ACl/GXMsiQ9f/58OgkAMHSAUueWW26xf/mXf3GuA1etWtvf/M3f2Jw5c+z555+ngwAAQwcAAMDQAQAAAEMHAAAADB0AAAAwdAAAAAwdAAAAMHQAAADA0AEAAABDBwAAwNABAAAAQ4fc6evrS9lftGiRPfDAA0XZ7rzzzqJpox9uffoGffpm8vWXLVuGoQeJ//3f/03Z//a3v120tn71q18V9d+Cfnj16Rv06ZvJ15epY+gBNnRFa8Xi//7v/4r6b0E/vPr0Dfr0zeTrY+gBN/R8BxAAAMIJhk6GTqSMPp8dfeYOGTpMtaFTQ0e/FPXpG/Tpm8nXD7ShDwwM2MyZMzO+p7Oz0znzr7y83GbMmGEPPvigDQ0NOceuXr1q69atc15bs2aNs/3hD3/IesyNV5182iBDR58sC33mDvqhzNDLysoSWybq6uqst7fXMc+xsTHbvn27LV682DnW2tpqHR0dife2t7dbc3Nz1mNuvOrk00Y6Q6eGDgAAgc/Q48aeL9OnT3d+1tTU2PDwcOJ1Ze5VVVVZj7nxqpNPG2To6JNloc/cQT/UNfR8Df3QoUNO1i4qKipSjimLj7+W6Zgbrzr5tJHO0Kmho1+K+vQN+vTN5OtHytD379/v3FktXkOf6G/jr2U6lsvruejk04buECcz17J8MtXV1c4Ax6M2/fRrP76hj36++wcOHPD986IfDf2urq6izPco6EfG0FtaWqy+vj7lpLMwZOjU0AEAIDKGfvjwYefEODe1tbU2MjKS2FfmXllZmfWYXzr5tJHO0Kmho1+K+vQN+vTN5OuH0tBlit3d3c7vBw8eTHspWFtbW8oStn5vamrKeszdhledbG3kYui1X/ta0fqWWhb69A36zJ3g6IfisrXky9cGBweda85HR0fTvi/+Xq/XiLvbmMrr0FdmuQ6fSBZ9siz0mTtk6IHI0N309PRYY2Nj4NvI1dAfmDXLAAAAQmfoDQ0N1t/fH/g2cjX07/7XfxHJok+WhT5zB33u5R403Ia+/N//vWhtUctCn75Bn7kTHH0MPeCGXv/VrxLJok+WhT5zB30MPeiG/kARDR0AAIIDhh5wQ1/9r/9KJIs+WRb6zB30MfSgG/p9//zPRWuLWhb69A36zJ3g6GPoATf0+7/yFSJZ9Mmy0GfuoI+hB93QH/jyl3WHGjoGACDiYOhBz9D/6Z/Mrt2xjkgWfbIs9Jk7ZOgYekAN3amhf/RRUdqiloU+fYM+cyc4+hh6wA3dOcv99GkiWfTJstBn7pChY+hBNnTnOvTf/Y6OAQCIOBh6wA29/j/+w+y3vyWSRZ8sC33mDhl6cA19YGDAZub4+NCJ3uvXo02n8vGpzr3cX3+9KP1LLQt9+gZ95k5w9EPxPHSv721tbbWOjo7Efnt7uzU3N2c95sarTj5tpDN052lrrteIZNEny0KfuUOGHqgMPW7WXt9bU1Njw8PDif2hoSGrqqrKesyNV5182khn6A987WtmL79sAAAQbSJt6BUVFSn7WgKPv5bpmBuvOvm0kc7QV95wg1lXF5Es+mRZ6DN3yNCja+gT/W38tUzHcnk9F5182ujr63PMXMvyydROm2bvP/JIYpD1UzUWP/bdP9FHP9f9rmtBJvro57u/ZcuWosz3KOiToQc8Q7/3xhvNfvpTIln0ybLQZ+6QoUfX0Gtra21kZCSxrxp2ZWVl1mNuvOrk00Y6Q39g/nyzxx83AACINqE0dJlid3d31ve2tbWlLGHr96ampqzH3G141cnWRi6Gfo9OovvhD4lk0SfLQp+5Q4YeTENPvhQt+ZK0wcFBKy8vt9GkB5ake6/Xa8TdbUzldeiLv/lNs7Vri9LHXA+KPn2DPnMnOPqhu1NcT0+PNTY2Br6NXA39bi25r1pFJIs+WRb6zB0y9HAZekNDg/X39we+jVwN/YG6OrNlywwAAKIN93IPGG5Dr1u82EwbkSz6ZFnoM3fI0DH04Br6f8+dazZvXlHaopaFPn2DPnMnOPoYesAN/c477zSbPZtIFn2yLPSZO2ToGHqQDd0ZQGXply7ROQAAEQZDD0OGvmiRWdINaohk0SfLQp+5Q4aOoQfM0L/97W+bLV1qduaM721Ry0KfvkGfuRMcfQw9DBn6ihVmJ04QyaJPloU+c4cMHUMPqqE7A7hmjdnRo3QOAECEwdDDkKE3NJi98QaRLPpkWegzd8jQMfSgGrpTQ1+3zmzfPt/bopaFPn2DPnMnOPoYehgy9CefNNu9m0gWfbIs9Jk7ZOgYelAN3RnA1laznTvpHACACIOhhyFDf/ZZPUydSBZ9siz0mTtk6ME09IGBAZs5c2bG90zGs8qn8nnoTg39xRfNWlp8719qWejTN+gzd4KjH1hDLysrS2yZaG1ttY6OjsR++3gm29zcXNCxqWgjY4a+d6/Zxo1EsuiTZaHP3CFDD2aGHjf2TNTU1Njw8HBif2hoyKqqqgo6NhVtpDN0ZwAPHDB77DEDAIDoEnpDr6ioSNnXMnf8Na/HpqKNvr4+x8zbXbXy6upqO7Jtm43W1yeiNy3JxKO4Qvbjm1966EdH/4CCTJ8/L/rR0O/q6irKfI+CfugNfaLj8de8HpuKNtJl6E4N/fhxs2uG7ifUstCnb9Bn7gRHnww9IBl6OkN3augffGB2112+9y21LPTpG/SZO8HRD72h19bW2kjSo0VVp66srCzo2FS0kc7QnQE8d87sttsMAACiSygNXabY3d3t/N7W1pZSd9bvTU1NBR2brDZyztA//dRs/nwiWfTJstBn7pChB8/Qky9bS758bXBw0MrLy210dNTZL8Y14pPRRq6G7tTQx8bM5szxvY+pZaFP36DP3AmOfujuFNfT02ONjY2BbyOvDF3ccovZ5ctEsuiTZaHP3CFDD4ehNzQ0WH9/f+DbyNXQEwO4cKHZxx8bAABEE+7lHjDSZuh33GF29iyRLPpkWegzd8jQMfQgGrpTQxfLl5udPOlrW9Sy0Kdv0GfuBEcfQw9Lhn7//WbHjhHJok+WhT5zhwwdQw+ioScG8OGHzd58kw4CAIgoGHpYMnQ9pe2VV4hk0SfLQp+5Q4aOoQfR0BM19M2bzfbs8bUtalno0zfoM3eCo4+hhyVD37bN7LnniGTRJ8tCn7lDho6hB9HQEwPY0WH285/TQQAAEQVDD0uG/sILZk8/TSSLPlkW+swdMnQMPYiGnqih9/SYbdrka1vUstCnb9Bn7gRHH0MPS4b+6qtmP/oRkSz6ZFnoM3fI0DH0IBp6YgAPHTJ76CE6CAAgooTe0MfGxmz9+vU2Y8YMmzlzpm3cuDFxzK9Hm07l41MTGboeFnPffUSy6JNloc/cIUMPp6Fv2rRpPHk95Px+/vx527Bhg23dutXZb21ttQ6dHX6N9vZ2a9YNWrIcc+NVJ5820hl6ooY+MGB2zz2+9h21LPTpG/SZO8HRD72hKytPRhnwrbfe6vxeU1Njw8PDiWNDQ0NWVVWV9Zgbrzr5tJE1Qx8cNLv9diJZ9Mmy0GfukKGH09CnTZvmLG0nU15e7vysqKhIeV3vi7+W6Zgbrzr5tJHO0BMDeOGC2YIFBgAA0ST0hr5y5UpnWVu1dBnmG2+8kTD0srKy694ffy3TsVxez0Unnzb6+vocM9eyfDLV1dXOEszAe+/ZeHrvRG/aj0dxhezHN7/00I+O/oEDB3z/vOhHQ7+rq6so8z0K+qE3dC2xr1692jFxnRjX1NRkc+fODU2GnqihCy3XX7niW99Ry0KfvkGfuRMc/chdtnbw4EHnzHJRW1trIyMjiWOqYVdWVmY95sarTj5tpDP0RA1daMldS+8+QS0LffoGfeZOcPQjZehHjx4d97wFjnGKtra2lCVs/a4MPtsxIePt7u4uSCdbG7kYesoA6qQ4nRwHAACRIxKGruV2bXV1dXby5MnE616vER8cN03pjY6OFqTj63XoQpet6fI1Iln0ybLQZ+6QoYc5Q/eLnp4ea2xsnJK2M9bQdWMZ3WDGJ6hloU/foM/cCY4+hu6BhoaGcd/sn5K2M2bouvXrtZvoEMmiT5aFPnOHDB1DL2Ey1tD1cJbXXqOTAAAiCIYecENPydCfeCL2GFUiWfTJstBn7pChY+jBMvSUGvrTT5u98IJvbVHLQp++QZ+5Exx9DD1MGfrPf272zDNEsuiTZaHP3CFDx9CDZugpA/jcc2bbttFJAAARBEMPU4a+Z4/Z5s1EsuiTZaHP3CFDx9CDZugpNfRXXjHL8jz1fKCWhT59gz5zJzj6GHqYMvTeXrOHHyaSRZ8sC33mDhk6hh40Q08ZwGPHzO6/n04CAIggGHqYMnTdp/7ee4lk0SfLQp+5Q4aOoQfN0FNq6B9+aHbHHb61RS0LffoGfeZOcPQx9DBl6B9/bLZwIZEs+mRZ6DN3yNDDZ+h69vnq1autoqLC2ZYvX25nzpxxjvn1aNOpfHxqygBevmz2rW8R9QAARJDQG/qCBQusu7vbMU9tO3bsSCxTt7a2WkdHR+K97e3t1nztsq9Mx9x41cmnjZwydDFnjtnYGJEs+mRZ6DN3yNDDZejl5eWOkSczffp052dNTY0NDw+nZPNVVVVZj7nxqpNPG+kMPaWGLubPN/v0U1/6jloW+vQN+syd4OiH3tA7Oztt69attmvXLhscHLS9e/c6zzMXWoJPRsYffy3TMTdedfJpI+cM/bbbFBkQyaJPloU+c4cMPVyGLhNfsmSJY+I333yzU0+/cOGCc6ysrOy698dfy3Qsl9dz0cmnjb6+PsfMtSzvNnRFbPFB/mz83/r6jh2Jff1MPs4+++yzz34490Nv6LW1tXbixInEvmrWy5YtC2+GXl9vdvw4kSz6ZFnoM3fI0MNl6PF6+USvyexHRkYSr6uGXVlZmfXYREGDF5182khn6NfV0L//fbPDh33pO2pZ6NM36DN3gqMfekPXZWpado9nwC0tLeNJbL2z39bWlrKErd+bmpqyHhMyXp09X4hOtjY8ZeiNjWYHDhDJok+WhT5zhww9XIau67p1fbfOdtemWnq8hu71GnEFCNIaHR0tSMf369DFxo1me/caAABEC+4U54Genp7xRLhxStrOmqG3tJi9+CKRLPpkWegzd8jQMfRsKMvv7+8vCUO/rob+s5+ZPfusL21Ry0KfvkGfuRMcfQw9YGTN0Hfu1C3oiGTRJ8tCn7lDho6hB8nQrxvA3bvNnnySjgIAiBgYetgy9H37zNatI5JFnywLfeYOGTqGHiRDv66G/sYbKvL70ha1LPTpG/SZO8HRx9DDlqEfPWq2Zg2RLPpkWegzd8jQMfQgGfp1A6jb3K5YQUcBAEQMDD1sGfqZM2ZLlxLJok+WhT5zhwwdQw+SoV9XQ9e94Rct8qUtalno0zfoM3eCo4+hhy1Dv3TJbO5cIln0ybLQZ+6QoWPoQTL06wbw6lWz2bNjPwEAIDJg6GHL0MW8eWYXLxLJok+WhT5zhwwdQw+KoV9XQxeLF5sNDxfcFrUs9Okb9Jk7wdHH0MOYoS9bZnb6NJEs+mRZ6DN3yNDDY+hlZWUTbsKvZ5WX1PPQxXe/a/buu0Q/AAARInIZ+vDwsM2fP9/5vbW11To6OhLH2tvbrbm5OesxN1518mkjrwx97Vqz3/6WSBZ9siz0mTtk6OE19KamJtuzZ4/ze01NjWPwcYaGhqyqqirrMTdedfJpI52hT1hD/+EPzX7zm4L7iloW+vQN+syd4OhHytDPnDlj1dXVif2KioqU41oCj7+W6Zgbrzr5tJFXhv7443ojkSz6ZFnoM3fI0MNp6I888ojt378/sR+vpScTfy3TsVxez0Unnzb6+vocM9eyvNvQFbHFB1k/T69dax/9z/8k9t3H2WefffbZD99+ZAx9YGDAFutyLh8y65LP0NvazDo7iWTRJ8tCn7lDhh4+Q1+9erUdOnQo5bXa2lob0b3Pr6EadmVlZdZjbrzq5NNGOkOfsIYuM7+WoRcCtSz06Rv0mTvB0Y+EoZ84ccKWTvAEsrbxTDZ5CVu/66S5bMeEjLe7u7sgnWxteM7QX37Z7Cc/IZJFnywLfeYOGXq4DL2urs6OHTt23eterxEfHBy08vJyGx0dLUinaNeh6z06MQ4AACIDd4rzQE9PjzU2Nk5J2zll6Lpk7dFHiWTRJ8tCn7lDho6hZ6KhocH6+/tLwtAnrKHrpjK6uUyBUMtCn75Bn7kTHH0MPWDklKH/7nex278SyaJPloU+c4cMHUMPhqFPOIB6MIse0AIAAJEBQw9jhv7RR2bf+Q6RLPpkWegzd8jQMfSgGPqENXSdfT9vXsFtUctCn75Bn7kTHH0MPYwZ+tWrZrNnx34SyaJPloU+c4cMHUMvfUNPO4C33mr2xz/SYQAAEQFDD2OGLmpqzJJuK0skiz5ZFvrMHTJ03wy904cHhmDoOdTQhW51+/vfF9QWtSz06Rv0mTvB0Z9UQ9ejQW+66aastzcFHzL0lSvN3nuPSBZ9siz0mTtk6P4b+g033OCYurYNGzbgzj4YetoBXLPG7K236DAAgIgw6TV03Qddz/yWqc+YMWPcczCdomToupf7668TyaJPloU+c4cMvTiGHkf19GnTpiUy9vim1yB3Q09bQ1+3zmzfvoLaopaFPn2DPnMnOPpTZuh69neyice36dOn+94BJ0+etBUrVjgrA9ri+PVo06l8fGraDP2pp2LPRSeSRZ8sC33mDhl6MQx9z549iSX3mTNn2vHjx4v6jz979qzNnz/fjhw5ct2x1tZW6+joSAkympubsx7zSyefNtIZetoBHNe2nTtZ0gAAiAiTauizZs1KZOVPKYOcBGSQr7766oTHampqbHh4OLE/NDRkVVVVWY/5pZNPG3ln6Dt2mP3sZ0Sy6JNloc/cIUP339Bl5DKs8+fPT1rEovZk6joBTysD9fX1iWXt5OV3oSXw+GuZjrnxqpNPG+kMPW0Nfdcus5aWgvqOWhb69A36zJ3g6E+qob9cYE3XC6rL9/b2OmY5NjZmW7dutZW6RvtagDFR0JHtWC6v56KTTxt9fX2OmWtZPpnq6mpngONRm35qf+iZZ8w2bEjsu4/nsh/fvP49+tHVP3DggO+fF/1o6Hd1dRVlvkdBP/S3fnWfNS9jLy8vD02GnnYADx40a2w0AACIBqE3dC1JJ9ephW5wI2pra20k6X7nqmFXVlZmPebGq04+baQz9LQ19MOHzR58sKC+o5aFPn2DPnMnOPqhN3QtUa9fvz5RNz906JBt3LjR+b2trS1lCVu/NzU1ZT0mZLzd3d0F6WRrIxdDT1tDf/dds1WrCuo7alno0zfoM3eCox+Jp61t2rTJOSlOS+0y9ytXrjive71GfHBw0NEaHR0tSKeo16F/8IHZXXcRyaJPloU+c4cMPTyG7je6fW3jFNWnc66hDw2ZLV7MYAEARAQM3QMNDQ3W399fEoaeNkP/9FOz+fOJZNEny0KfuUOGjqGXIjnX0D//3GzOnILaopaFPn2DPnMnOPoYesANPW2GLr71LbNLl4hk0SfLQp+5Q4aOoZe6oWccwIULzT7+mE4DAIgAGHqYM/QlS8w+/JBIFn2yLPSZO2ToGHqpG3raGrq49149O9ZzW9Sy0Kdv0GfuBEcfQw9zhq7BffttIln0ybLQZ+6QoWPopW7oGQfwkUfMenvpNACACIChhzlDb24227+fSBZ9siz0mTtk6Bh6qRt6xhr6k0+a/eIXntuiloU+fYM+cyc4+hh6mDP0bdvMnnuOSBZ9siz0mTtk6Bh6qRt6xgF85hk9wo1OAwCIABh6mDP0F14we/ppIln0ybLQZ+6QoWPopW7oGWvov/yl2RNPeG6LWhb69A36zJ3g6Ife0MvKyibchF/PKi/J56GL114z+9GPiGTRJ8tCn7lDhh4OQ09Ha2urdXR0JPbb29utWZd6ZTnml04+baQz9IwD2Ndn9oMfsKwBABABIm3oNTU1Njw8nNgfGhqyqqqqrMf80smnDU8Z+jvvmN13H5Es+mRZ6DN3yNDDYejTp093tptvvtk2b95sly9fdo5VVFSkvFdL4PHXMh1z41Unnzb6xrNtmXm766z12bNnOzWV+CDrZ2J/YMAuL1mS/niWfffPfP8e/ejqd3V1+f550Y+G/pYtW4oy36OgH6mT4pQNy9DXrl2bNnuPv5bpWC6v56KTTxueMvRz58xqa4lk0SfLQp+5Q4YeLkOPZ8Hl5eWBy9DTGXrGAbxwwWzBAgMAgPATOUO/ePGiff3rX3d+rx3PXkdGRhLHVMOurKzMesyNV5182vCUoV+5YjZnDpEs+mRZ6DN3yNCDb+h1dXXW29vrZL8yc10m9tRTTznH2traUmrS+r2pqSnrMSHj7e7uLkgnWxu5GHrG69DFLbeYffaZp77jelD06Rv0mTvB0Q+9oT///PO2bNkymzZtmt14443W0tKSOOb1GvHBwUFn2X50dLQgnaJfhy605H7+PJEs+mRZ6DN3yNCDbejFoKenxxobG6ek7bxq6OKOO8zOnmXQAABCDobugYaGBuvv7y8JQ8+aod9zj3P5GpEs+mRZ6DN3yNAx9BIi7xr66tVmx455aotaFvr0DfrMneDoY+gBN/SsGfpDD5kdOkQkiz5ZFvrMHTJ0DL2UDT3rAOqs+VdfpeMAAEIOhh72DH3TJp3FRySLPlkW+swdMnQMvZQNPWsNfetWXbvnqS1qWejTN+gzd4Kjj6GHPUPX41mTHtFKJIs+WRb6zB0ydAy9BA096wA+91wsSwcAgFCDoYc9Q1f9XHV0Iln0ybLQZ+6QoWPopWvoWWvor7xi1tzsqS1qWejTN+gzd4Kjj6GHPUN/883YtehEsuiTZaHP3CFDx9BL19CzDqDuEnf//XQcAEDIwdDDnqErmlu+nEgWfbIs9Jk7ZOjhMfTdu3dbWVlZYt+vR5tO5eNTs9bQ9aQ1PXHNA9Sy0Kdv0GfuBEc/MoZ++PBhW7VqVYqht7a2WkfSNdrt7e3WfO0EskzH3HjVyacNzxm6noVeXU0kiz5ZFvrMHTL04Bv6qVOnbNmyZTY6Oppi6DU1NTY8PJzYHxoasqqqqqzH3HjVyaeNdIaedQA/+8zsllsMAADCTegNXUvYtbW1NjIy4uwnG3pFRUXKe7UEHn8t0zE3XnXyacNzhi7mzDEbGyOSRZ8sC33mDhl6cA397rvvtpMnTyb2kw09+Xf3a5mO5fJ6Ljr5tNHX1+eYuZblk5k9e7ZTU4kPsn669z+bO9fef/vttMfT7bt/5vv36EdXv6ury/fPi3409Lds2VKU+R4F/dAbugxyoi1SGXptrdm5c0Sy6JNloc/cIUMPrqFnyoCTl+KFatiVlZVZj13vl9508mkjnaHnNIB33WX2/vsGAADhJdKG3tbWlrKErd+bmpqyHhMy3u7u7oJ0srXhW4Z+331m77xDJIs+WRb6zB0y9HAautdrxAcHB628vNw5a74QnUm5Dl18//u6bi/vvuJ6UPTpG/SZO8HR505xHujp6bHGxsYpadtThv7YY2avvUYkiz5ZFvrMHTJ0DD2ZhoYG6+/vLwlDz2kAN24027uXgQMACDEYesDwlKFv2WL24otEsuiTZaHP3CFDx9BL1dBzqqHrxLtnn827LWpZ6NM36DN3gqOPoUchQ9fZ+Nu3E8miT5aFPnOHDB1DL1VDz2kAf/ELsyefpPMAAEIMhh6FDH3fPrMf/5hIFn2yLPSZO2ToGHqpGnpONfTeXrNHHsm7LWpZ6NM36DN3gqOPoUchQz961GzNGiJZ9Mmy0GfukKFj6KVq6DkN4IkTZvfeS+cBAIQYDD0KGfqZM2ZLlxLJok+WhT5zhwwdQy9VQ8+phq77wy9cmHdb1LLQp2/QZ+4ERx9Dj0KGfumS2dy5RLLok2Whz9whQ8fQS9XQcx7A2bPNPv+cDgQACCkYehQydKGl+U8/JZJFnywLfeYOGXowDX3Hjh22dOlS5/nlFRUVtnLlSjt16pRzzK9nlZf889DF4sVmw8N59R21LPTpG/SZO8HRD72h19XVWW9vr2Oe2jo7O+3WW291jrW2tlpHR0five3t7dbc3Jz1mBuvOvm0UXCGPt4Pdi2QIZJFnywLfeYOGXooltynT5/u/KypqRlPWv+UtQ4NDVlVVVXWY2686uTTRjpDz3kAV60ye/ddAwCAcBIpQx8bG7OdO3c6S/BCS/DJKIOPv5bpmBuvOvm00dfX55i5svhkqqurnSWYeNSmnxPtj8rQjxxJe3yi/fiW6/vz3Uc/vPoHDhzw/fOiHw39rq6uosz3KOhHxtDLysqcTXXq8+fPJ16b6H3ZjuXyei46+bSRLkPPuYbe2Gh28GBefUYtC336Bn3mTnD0I5WhX7lyxTlJbsmSJYHL0NMZes419Mcf10zIq7+oZaFP36DP3AmOfiRr6DrjXdTW1trIyEjiddWwKysrsx5z41UnnzbSGXrOA9jSYrZrlwEAQDgJvaEvX77cjhw54mS/ly5dcmro9fX1zrG2traUmrR+b2pqynpMyHi7u7sL0snWhq8Z+nhb1tlJJIs+WRb6zB0y9GAa+u7du51L13Rm+6xZs5zLwi5evOgc83qN+ODgoJPlj46OFqQzqdehd3XpOrm8+o5aFvr0DfrMneDoc6c4D/T09FijTjKbAjxn6C+/bPaTnxDJok+WhT5zhwwdQ4/T0NBg/f39JWHoOQ/gr39ttn49gwcAEFIw9IDhOUN//XWzRx8lkkWfLAt95g4ZOoZeioaecw39rbfMvve9vNqiloU+fYM+cyc4+hh6VDL0994zW7mSSBZ9siz0mTtk6Bh6KRp6zgN4+rTcnw4EAAgpGHpUMnTdwKamhkgWfbIs9Jk7ZOgYeikaes41dF0zP29eXm1Ry0KfvkGfuRMcfQw9Khn61atms2fHfhLJok+WhT5zhwwdSy0tQ89rAJWhX7u7HQAAhAsMPSoZulANPelhMETK6PPZ0WfukKFDiRh6zjX0mPvHznbPEWpZ6NM36DN3gqOPoUcpQ9d16LoenUgWfbIs9Jk7ZOhYamkZel4DqDvF6Y5xAAAQOkJv6J2dnbZs2TLncaczZsxwHlU6NDTkHPPr0aZT+fjUvDJ03ctd93QnkkWfLAt95g4ZetAMXc9C7+3tdcxzbGzMtm/fbosXL3aOtba2WkdHR+K97e3tzvPSsx1z41UnnzbSGXpeNfTx4MF56lqOUMtCn75Bn7kTHP1ILrlPnz7d+VlTU2PDw8OJ15W5V1VVZT3mxqtOPm34kqHfe6/ZmjVEsuiTZaHP3CFDD76hHzp0yMnaRUVFRcoxZfHx1zIdc+NVJ5820hl6XgOo5fwVK8y2bTMAAAgXkTL0/fv326JFixI19LKysuveE38t07FcXs9FJ582+vr6HDPXsnwy1dXVzhJMPGrTz0z77x87Zhduu80+3rAh6/vjWz76+eyjH179AwcO+P550Y+GfldXV1HmexT0I2PoLS0tVl9fn3LSWRgy9Lxq6HE++cTsnntUtM/4NmpZ6NM36DN3gqMfCUM/fPiwc2Kcm9raWhtJunOaMvfKysqsx/zSyaeNdIZ+p9dHol64YHbXXWZJJ+W5oZaFPn2DPnMnOPqhN/SDBw+mvRSsra0tZQlbvzc1NWU9JmS83d3dBelkayMXQ893AFP4+GNdBmC2Y4cBAECwCb2hqyY90Sa8XiM+ODjoXNc+eu1BJ4G5Dn0i1J40du4kUkafz44+c4cMPVr09PRYY2PjlLTtSw3djZb9ly41e/75lJepZaFP36DP3AmOPobugYaGBuvv7y8JQy84Q4+j6+HvuMNs1y4iZfTpG/SZO2ToMNmGXlAN3Y0u57v9drOXX6ajAQACBoYecEP3LUOPc+6cTr8327OHSBl9+gZ95g4ZOkyWoftSQ3dz9qzZbbfZO088UdR/C7Wy8OrTN+jTN5Ovj6GToU/Mhx/a59/6lq6tI1JGn75Bn7lDhg7FNnRfa+huDh82W7jQ7N136XgAgBIHQydDzxwJHjqkx8LFautEyujTN+gzd8jQoTiGXpQa+jUStZqXXordJvbajXR81y/250d/0vXpG/Tpm8nXx9DJ0HOLBH/6U7Pvf9/s88+JlNGnb9Dns5Ohg9+GXtQaejIy8oceMvvJTxgEAIASBEMnQ889EtSSu5bek+4mR6SMPn2DPp+dDB18MPRJqaEno5PjdJLcm28WR7/Ynx/9SdGnb9CnbyZfH0MnQ88/EtRlbNXVZgMDRMro0zfo89nJ0CeXgXHzmTlzZsprfj3adCofnzppNXQ3r70Wu0Vsls8LAACTQyQM3f0c9Ditra3W0dGR2G9vb7fm5uasx9x41cmnjZLK0OM8+6zZypVmly8TKaNP36BP35ChT66xJ1NTU2PDemzoNYaGhqyqqirrMTdedfJpI52hT3oN3c369WaPPablhuLoF/vzo89nR5+5ExL9SBt6RUVFyr6WwOOvZTrmxqtOPm2UZIYurlwxu/9+T/d8JxIny0IfffqGDN0XQ3fvJ7+W6Vgur+eik08bfX19jplrWd5t6IrY4oOsn5O9/8pLL9lnixcr2piS9tlnn3322SdDJ0P3K9L8/e/N5s41q68303kAugHNz35m9txzZnv3mr3+utmxY2bvv2/20Udmly7ZwIkTROJkWeijT9+QoRdu6LW1tTYyMpLYVw27srIy6zE3XnXyaSOdoU95DT0ZmbYe5vLKK2Yvvxw7aW7rVrONG80efTS2NK8b03znO2Z6NKv+rQ8/bPbee6Xx+dHns6PP3AmwfqQNva2tLWUJW783NTVlPSZkvN3d3QXpZGsjUBm6Bz54661YBq/L31avNjt4kHvFk2Whz9yhb8jQMxu5exNerxEfHBy08vJyG7329LFIXofuJzLxAwfM7rvP7PbbzV54wfcnuwEAhB3uFOeBnp4ea2xsnJK2w5ShT6iv5XfV4P/7v82eflrRE5E4WRb6zB30ydCLQ0NDg/X395eEoZdUDd1PfV2f39pqtmCBmYInnVAXpM8fcX36Bn36ZvL1MfSAEfoM3c2lS2a7d5stXRrL2n/wA7Mf/chs3brYyXZPPmnW0mK2bVvsWnjdea+z06nNf6rnt+u+82NjRPpkWegzd8jQMfTSNvRQ1NBzQXei00mIOpNe95Hfty92OdwvfhF7nKtOrpORy9Bl7Dq7/u67ze65x2zePLPvfc/s5z83++1vzf74RyYSAIQODJ0MPfyR7MWLscvpZPS6dO7WW81WrDDbsiV2Zv3HHxPp89nRZ+6QocPUGnpoa+jF1Nftat95x2znTrNHHokt5d9xh6KjWIb/y1+aHTlidvp0bMk/av3D3EGfuRNIfQydDJ1IWcv5J0/G6vFalldtfu3aWN1eN8CR4S9fHjN/3QFPgcC+fTasewh89hn9T5aFPp+dDB0KN/TI1NCnkgsXYoav+v1LL5lt365LHWJGr9vdavleRq+xOXPG05PnAAAKBUMnQydSLkRfGfrx47Gb4ej6ed0YZ/58M51hr2xftXsFBGRZ6KNP35ChQyZDp4Zegvo6ya63N/ZwGi3da4xuuy1Wo9eJeFqy1zgePhx7WM358xNm9dRB0Y+iPn3jXR9DJ0MnUi62vsxaWbxq7srkdcOc9evNHnww9rAa3Tynqir20JqVK2O1+sces4s6I1+X6u3ZE3vgjbJ9ncz3wQexG+/o7P0ClvfJstBn7pChQwkZOjX0kKCb38ikdevbN94w27EjtoSvev3mzWZ6aM9DD8Xud68gQM+f1/X1s2fHfioo0Pbd78auuVdQoL/ftCl2C10FE1oZ0FPwNId0sx6dza8AYWAgdotdrRQU8SQ/ACguGDoZOpFykPWVoetBNrrETrfH/d3vzPQUOwUF+/frwQOxVQE9ylYrAzp57/HHnefWX9alegoQdAMe1f4VENxyi9mcObHzABQ0LFsWCxL0Uyf//fCHseBCKwwZ7tT3id6nfd3859VXYyUI3dRHd+5TAPHhh2Z6dLBWGTzcyY+5Q4aOPoYeOkOnho6+7/oy2E8/ja0YnDoVCxL0XmX0v/lNzKB//euMd+ob1mV+ukWvTF8BgFYMdD7BqlWxu/cpmKipiZ1foACisjJ2xYCCCr2uR+oqWNV7VYZQ4LFmTezWv48+ahdVntCKhdqYaHvssT9tehaA2tOVCT/+8Z+CEa1eKCBRkKNzGxSUKOgZ//x/UACjcyAUCOnfpdUN/RsVHOnfq6sdtMqhcoiCpieeiPWJgij1j54eqL5SYPXmm7HzJRTQHD3qbMfVvoIbrcicOKH/pWPnU6i/f//72NUSZ8+anTsXWz3RMZ1cqbsc6j4KWUot1NCjqY+hk6ETKaM/9dqXL8eCCJmWMneZmIxN2bxMT6UBmWFfn3P54McyYwUYug3wRJsMNb7pboAqOcQNNx6MyIgVkKgM8eKLMcOWcY8b+Gh9fexvtGmlQUYvw1f5QisSP/2p2VNPxcohMmd9DxUoKGBQAKNAQqsUCiIefjh21YMCGgUl49qfq0yi4EbBilY+FABppaSuLqaleyAoCFFgs2hRbOVEwY7ucqjzLRQAKRBSEKTVFL1HJ14uWeL83ecKlKSrNhQM6T96lWL0ORQI6XMp0FHwo3M59Dn0+bV6o23Dhti/S5uCFfV3fBsPzkalq37QLZavBUFOMKcAKB78KPBR3yow0hUf6nP9/6X+162bNRY6NyQ+Zhqn8SDoI5WZNLY6ZyQ5ENIKlOZBPBDSpaQKdHRM+1qligdD8U0rQQqMtGlOjW9n1K5eV8CquaZHVutEVpWcNP80D7VypJUv3VhK+5qfCnRzOGeFDD2i8Dx0APCMDCZuOPp/Y2goZlQyN5U4lPnL/GSCb78dK8XIIGWUWjmQgcpItdoQD3j0/4s2ZYEyYG0yVwU/8U3mLCOXYT//fCIIcgxdOvHgR4av1Q+tqihoUHCgYEGBgx6upOBH53nEV1UUYCgIUnCjMo8Cj+RASCeJKjiJB0L33hsLghToxMtD8WAovl0LcJxNZSVtCooWLoz9jVaDFAxpv7o6pqUASQGRgi4FTPHgKb6SdPPNf3pNgZbeI039jf5W57XoHhXx81r0uzYdc2/6G20+PUAKQ59CWscnfoe+BNdoH/8yNGuCk6Gjb9RB0WfulKy+MnWZsE4iVVClUoiy+vHg6rSCqXi2r+2TT2KbAi/3pr/RRoYefGrGo8NhLftcY2g8wq7Scloehk4NHf1S1Kdv0KdvJl8fQ59CKioqXEHf1eteI0NHnywLfeYO+mToJU5ZWVlOr4m+vj7HzJ9XzSqJe+65xxlENjY2NrZob8t0HgGGHpwMHQAAwA8wdB+pra21EV02cQ3V0Ct19iQAAACGHhza2tqcM9vj6PcmXcMKAACAoQcHL9ehu+nu7nZq68XYFGAUSxv9cOvTN+jTN5Ovf1j3GcDQg4v7rPegaKMfbn36Bn36pvT1MfQSI9+IrFS00Q+3Pn2DPn1T+voYOgAAQAjA0AEAADB0AAAAwNABAAAAQw8Lflzulo7Ozk7n9oHl5eU2Y8YMpw3d8KYY7N69O+2tbgvh5MmTtmLFCueue37feU99sXr16oT28uXL7Yyeu1wgAwMDNnPmzKKN80T6fo71RPp+jnU6fb/GeiL9Qsc6U//6MbaZ9P0Y21w1vI5tNv1CxzaTvh/f4x07dtjSpUsdfWmsXLnSTp065cv4ZtL283uLoZcAXh67mit1dXXW29vrTMixsTHbvn27LdZzg31GZ2quWrXKd0M/e/aszZ8/344cOVKUvl+wYIFz7b/6R5u+eIU+8U59EN+KMc7p9P0a63T6fo11On2/xjqdfqFjnal//RjbTPp+jG0uGoWMbSZ9P8Y2k74f3+NkfW0y2lv1XHQfxjeTtp//R2PoJYCXx64WwvTp033VU6SpCHN0dNR3Q9eX5tVXXy1aXygq1hepGP3j7gu/xzmXvi7k3zKRvp9j7f57v8farV+MsY7/fbG+w5k+nx/zNFmjGN/juH6xvsdx/WJ9j4s5vsUYWwy9BJjMh7ocOnTIiQj9QstOyfew99vQ9aXRfwZailKf1NfX+1aOiC93bd261Xbt2mWDg4O2d+9ea2hoKIqh+D3O2fq60LF26/s91u6/93us3fp+j3Vy/xbjO5xp/Pz4HidrFON7nKxfjO9xsr7fY6tMeefOnc4yud/j69b2c2wx9BIgn8euFsL+/ftt0aJFvtbQ7777bqc2VqzPPW3atJTlKH1pVX/yC335lyxZ4nz5b775ZqcOd+HChaKMq9/jnOlv/Rhrt77fY+3+e7/H2q3v51i7+9fvsc00fn6MrVvD77F16/s9tm59P8c2Xq5Rnfz8+fO+ju9E2n6OLYYekQy9paXF9+w2eYK6Nz8N3d03Wl7zC2UlJ06cSOyrTpbvM4hLLUP3a6wnCkj8HOuJDN3PsXbr+zXWE/Wvn2Obafz8GNuJNPwc24n0/RzbifT9/h5fuXLFqcMrSPB7fN3afo4thl4CFPuxqzrRRdHxVK02FIJObEmuXYkbbrjBN/2JalXFqqH7Pc4T9bWfY51tLP3O0P0ea7e+H2Odrn/9GttM4+fH2Oaq4XVs0+n7Nbbp9Iv1PY4HHcX4Pzo5oPHre4uhlwDFfOzqwYMHfc/KJ9PQ1Rfr169P/BtUX9q4caNv+rq8Rct18ag7HiUXoy/8Hme3vt9jPdmG7vdYu/ULHetM/evH2GbS92Ns89HwMraZ9P0Y20z6fnyPpaGz8PX3ly5dcurccY1CxzeTtp/fWwy9BCjmdejFXhIvtqGLTZs2OSfTKKLVfwpasvIL9bP6W9raVIMrtIaerr/9Gud0+n6Nda4afly25tb3Y6zT6Rc61pn614+xzaTvx9jmo+HXvPFzbDPp+/E91vX3OhlNmf2sWbOck/guXrzoy/hm0vbz/2gMHQAAIARg6AAAABg6AAAAYOgAAACAoQMAAACGDgAAgKEDAAAAhg4AAAAYOgAAAGDoAFAS9PT02Ne//nXnzmEToTtq6bjeBwAYOgCUMLolqG51+cwzz6S8rvtn6/XNmzfTSQAYOgAEAd0bW+b92muvOft79+5NPDcaADB0AAgIelCHng+tZ0zrSVT6efvtt9vly5fpHAAMHQCCxPnz562qqsrJzOfMmTOpj/0FwNABAHxCBo6hA2DoABBg4kvuepa1To5jyR0AQweAAHL//fc7mfmvfvUrZ3///v3Ovl4HAAwdAAJA/LI1XaaWjDJ1va7jAIChA0AJE7+xjG4gMxG64Qw3lgHA0AEAADB0AAAAwNABAAAAQwcAAAAMHQAAAEMHAAAADB0AAAAwdAAAAMDQAQAAMHQAAADA0AEAAABDBwAAAAwdIOhfxi98IbAbAGDoAJBk6HxuAMDQATB0PjcAhg4AGCOGDoChAwCGDgAYOgBg6ACAoQNg6B4pKyvD0AEwdACIiqFneh+GDoChAwCGDgAYOgCGnqsZd3Z22pw5c6yiosLq6+vtwoULKUY9NjZm69evtxkzZjjbunXrnNfi70neMHQADB0AcjX0NWu8bxMYeltbm128eNGuXr3qmHtTU1OKoW/evNmeeeYZ57i2nTt32qZNm8jQATB0ACjI0I8e9b5NYOjJKPOeOXNmyrGbbropkZGLK1eu2I033oihA2DoAFCQofvIRGY8ffr0lGOZ3oOhA2DoAFCChq5M/Bvf+EbKMWXj7gw9/h4MHQBDB4ASMfSTJ086v6s+/uKLL9qGDRtSjFo19Pb29kQNXfX05Br6rFmz7NSpUxg6AIYOAFNp6HV1dVZeXu6cwa6z2S9fvpxi6MrOdWa73qNNvytLj7Nr167EGfAYOgCGDgBTZOhB/NwAgKEDYOhJKOPG0AEwdAAIuKHzuQEwdADAGDF0AAwdADB0AMDQASJo6EHdAABDBwAAAB/4f5ffaYggM7JSAAAAAElFTkSuQmCC",
      "text/plain": [
       "BufferedImage@4d357a79: type = 2 DirectColorModel: rmask=ff0000 gmask=ff00 bmask=ff amask=ff000000 IntegerInterleavedRaster: width = 500 height = 400 #Bands = 4 xOff = 0 yOff = 0 dataOffset[0] 0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(svals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which shrinks a little too fast. Lets look at it on a log-log plot instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGQCAYAAABYs5LGAAAsAklEQVR42u2dfWxV9f3HSQgxhBBCYhayLGZZYogxi3+YEP8wCyFbmFvIHsw25IcKAkWY0CpP+/Hg0CFWHc9KKetcdZVfxUFBqDAf2ITCsKIVHCpCERGRh4o8Fkop7x+fc3a73lLKvbf3nnvO97xeyU17Tm+/59vvPZ++v+/z+T50EwAAAESebjQBAAAAgg4AAAAIOgAAACDoAAAAgKADAAAg6ADgB1G3bim/Onq/y7zxxhuaMmWK7rrrLvXp00c33HCDevTo4X3/i1/8QqWlpWppaeEmAkDQARD0qLfNHXfcoTNnznAjASDoANESdNrm6tekSZO4kQAQdIBwiRb8l5tuukmzZ89WTU2NmpubvXN79uzRkCFDktrsW9/6Fo0FgKADREvQr/V+EzwTv+985ztenrlnz54aOHCg9/21yr9WWZ3Vqf3PLIc9d+5cffe731X37t1b3/fhhx9q9OjRnign6nPrrbdq5syZXX5EfurUqaQ6WPkAgKADOCHogwcPTvuxfTYE/c4777zqveXl5Z64X6sOJuwnTpzIuM2s89K2vP79+3MjASDoAOES9EwGxVVUVFwlsg0NDTp27FjOBb39q66uLul42rRpnotfsmTJVeczpaysLKmsxx9/nBsJAEEHiL6gt3fnNt0rE3HO5HdWrFiR9Ah9+PDhST9vamryzpuotz3/ve99L6P2shx63759k8phlDsAgg7ghKBbfrp9XjsoQW9PW7Ht7JVJ3nvv3r369re/3VqGXWv37t3cRAAIOkD4BD2T92cqwLkQ9M5y510Z0W8j3W00e9uR7fZ4HwAQdABnBL29Q8+noPfu3Tvp54npZl2hpKQkqaNw22236eDBg9w8AAg6gFuCPmDAgKRzlmdOpfz2U9rMBe/YsaNLgn733Xcn/fzFF1+86j3nzp3zptRdD5ueNnTo0KTyxo0b15qXBwAEHcApQbf5523P2fxvy6Nv2LCh0/JvueWWrE11S1BbW3vVoi+VlZWeU7c6bdq0yZu2lsrfanPbWUkPAEEHiI2g25xuW8DleoLXduEXw6Z7tX+PjRrviqAbNo3ONlLpqgCnI+YIOgCCDhB5QTcOHTqk++67T7169fIepf/gBz/Q5s2bk95rufb2PPzww955E+AHHnjA6xx0VdCNAwcOeDulWTrAyrfOhF3j+9//vgoKCjwnj6ADIOgAkAI2paut4KWStwYABB0A8oit4f7cc895q8MZH3/8sefS2wq6PQoHAEDQAcIcjNd5HD1ixAgaCQAQdICwM2zYMN1+++3ePHDLVdvLVlWzaV/r1q2jgQAAQQcAAEDQAQAAAEEHAAAABB0AAAAQdAAAAAQdAAAAEHQAAABA0AEAAABBBwAAQNABAAAAQY8TdXV1SceTJ0/W73//+5RfhYWFab0/rq8otVM+65rra2ez/GyUlWkZmfxeOr8zatQo4pa47vKrqKgIQQ+Sf/7zn0nH9iGkwz/+8Q8a0bF2ymddc33tbJafjbIyLSOT30vnd/7yl78QtMR1l0lXTxD0PAv6Z599RiM61k75rGuur53N8rNRVqZlZPJ76fzOv/71L4KWuEbQ4yboAAAACDoOnZ48dcWh49CJaxw6gh4GQSeH7l47kUMPrixy6MS1y3VF0HHo9ORx6Dh0HDqxgkOHoAUdAAAAQceh05Onrjh0HDpxjUNH0MMg6OTQ3WsncujBlUUOnbh2ua4IOg6dnjwOHYeOQydWcOgQtKADAAAg6Dh0evLUFYeOQyeucegIehgEnRy6e+1EDj24ssihE9cu1xVBx6HTk8eh49Bx6MQKDh2CFnQAAAAEHYdOT5664tBx6MQ1Dh1BD4Ogk0N3r53IoQdXFjl04trluiLoOHR68jh0HDoOnVjBoUPQgg4AAICg49DpyVNXHDoOnbjGoSPoYRB0cujutRM59ODKIodOXLtcVwQdh05PHoeOQ8ehEys4dAha0AEAABB0HDo9eeqKQ8ehE9c4dAQ9DIJODt29diKHHlxZ5NCJa5friqDj0OnJ49Bx6Dh0YgWHDkELOgAAAIKOQ6cnT11x6Dh04hqHjqCHQdDJobvXTuTQgyuLHDpx7XJdEXQcOj15HDoOHYdOrODQIWhBBwAAQNBx6PTkqSsOHYdOXOPQEfQwCDo5dPfaiRx6cGWRQyeuXa4rgo5DpyePQ8eh49CJFRw6BC3oAAAACHrUBf3YMZ0sKaER6cnj0HHoxDV1RdAjLehNTTo9bpy0bh0NeR3ItYXj2uTQU4McOnGNoMdN0K/wxQcfSAUF0q5dNCY9eRw6Dp24pq4IelQF3eOjj6TRo6WjR2lQAABA0KMo6K29vo0bpcmTpQsXaFR68jh0HDpxTV0R9KgJelJeZulSaf58GvV67URd83ZtcuipQQ6duEbQ4+zQjYsXpenTpdWraVh68jh0HDpxjUOPpqDX19erX79+rcdr165Vnz591KNHD/Xv3987zuf5XAn6VZw4IY0dK73/PtEOABBjIinoxcXFGjBggLp1+291Bg0apKqqKu/7iooKDRkyJK/nA3HoCfbs8QfJHT7MHU1PHoeOQyeucejREfSysjK1tLQkCXqvXr28c4Z97d27d17Pt8ece3l5uUpLS5POjxo1ysufJQLavnZ2/Ne//rXDn9cvXy4VFan27bfTKs/VY2unqNR3wYIFebv+te6nMJZvOcl81SeT+ymd+i5atCjW8RqW+zWbx9m4X9M9jnQOva2gd+/ePeln9gg8n+cDdegJ/vQn6emnpcuX6cnTk8eh49CJaxx6NAW9Z8+eHTrlfJ3PlaB3yqVL0qOPSitXEvkAADHDGUEfPHhw64C0devWebntfJ7Pi0M3Tp6UbHnY2lp68tQVh45DJ65x6NET9FWrVnn5bDvXt29fT1zzeT5Xgp7S3MZ9+yw5L33xRWwDn/mq4bg289BTg3noxHXsBT2KBLYf+ttvSxMnSmfP0pOnrjh0HDpxjUNH0LNFXV2dJ+aVlZVJ5wsLC72eXOLDt69ZOy4vV+OsWfpnrsrnmGOOOeY4NMcIuqsO3bDBeo89Jr30Ej156opDx6ET1zh0BD1Mgp52Xub0aem3v5W2bYtVO5NrC8e1yaGnBjl04hpBx6GnxoED/iA5+0pPnrri0HHoxDUOHUHPv6BnzNatvlM3xw4AAM6BoMfBoSeoqJAef9zPrdOTp644dBw6cY1Dh/TJ1ih32wgm41GQly+rceZMfTFnjvOjPhMb5kShvpY/zdf1u3Q/BVx++69B1ieT+ymd+paUlDgdj1G5X7N5nI37lVHuOPTOsXnpNj/d5qnTk8eh49Bx6MQ1Dh3yI+hZwVaQs0Fy9fV8IAAAjoCgx82hJ7C13m3Nd1v7nZ48Dh2HjkMnVnDoEKygZ3Vu48sv+7uz2S5tjsF81XBcm3noqcE8dOIaQcehdw3bN/2pp6SyMnryOHQcOg6dWMGhQ5CCnnUaG6WiIumtt/hwAAAiDIIeENmatmY9+WxPe3hnzRpdGjlS+vRTZ6a5JBxPFOprU3Hydf1c3E+5Kj/xykd9Mrmf0qlvdXV1pOPNlfs1m8fZuF+Ztua4Q89ZXub996WxY6UTJ5xoZ3Jt4bg2OfTUIIdOXOPQYyjoOc3LrF4tTZ8uNTdHvp3JtYXj2uTQU4McOnGNoMdQ0HPOvHnSf1atAgCA6ICg49CTuXBBmjRJ2riRnjwOHYeOQydWcOiQK0EPJC9z5Ig0erT00UeRbWdybeG4Njn01CCHTlwj6Dj03LFrl1RQIDU00JPHoePQcejENQ4dsi3ogfLqq9K0adLFi3xwAAAhB0EPiDDPQ+/s+OyTT+rIjBnMQ2ceet7LZx46x8xDZx66Uw498LxMU5M0daq0fn2k2plcWziuTQ49NcihE9c49BgKel5ySMeP+/l0y6tHBHJt4bg2OfTUIIdOXCPoMRT0vLF7tz/y/ehRPkQAgBCCoOPQU2fDBmnyZH+uOr1jHDoOHYdOXOPQEfTMP4C855CWLpXmzw99O5NrC8e1yaGnBjl04hpBx6EHj01hs/Xeq6roHePQceg4dOIah46gZ/oBhALbkc12ZrMd2gAAIBQg6Dj0zNizxx8kd/gwvWMcOg4dh05c49AR9Mjl0Nvy5ptSUZHU2Bi6dibXFo5rk0NPDXLoxDWCHiGiulLc9Y5Pz5un41On6rP9+0NRH1aKY6U4VopjpThWikPQI+HQQ8elS9KsWdLKlXy4AAB5BEGPmKCHMod08qQ0bpxUWxuaKpFrC8e1yaGnBjl04hpBj6GghzaHtG+fNGqUdOhQKKpDri0c1yaHnhrk0IlrBB2HHrY/Tpo4UTp7lp48Dh2HjkMnrhF0BD3SlJdLc+dKly/zYQMAIOgIemR7qDZI7rHHpBUr6Mnj0HHoOHTiGkFH0K9FJHJIp09Lv/2ttG1b3qpAri0c1yaHnhrk0IlrBB2HHl4OHPAHydlXevI4dBw6Dp24RtAR9Aizdav00EO+YwcAAAQdQY9mD9WjokJ6/HGppYWePA4dh45DJ64RdAQ9QZRySB422v2JJ6QXXgj0suTawnFtcuipQQ6duEbQI4Sra7mncrxl40ZdfPBBafNm1nJnLXfWcmctd9ZyZy13HHqkOXjQHyRXX89NAACAQ0fQI5dDb8s77/hrvtva7zmGXFs4rk0OPTXIoRPXCHoMBT1yOfT2vPyy9Oij/gI0OYRcWziuTQ49NcihE9cIOg49etgguaeeksrK6Mnj0HHoOHTiGkGPr6A7QWOjVFQkvfUWNwQAAIKOQ480hw9Lo0dLn35KTx6HjkPHoRPXCHr8BD3yOfS2vPeeZNPZTpzIetHk2sJxbXLoqUEOnbhG0HHo0Wf1amn6dKm5mZ48Dh2HDsQ1gh4fQXeSefOkZctoBwAABB2HHmkuXJAmTZI2bqQnj0PHoQNxjaDHQ9CdyqG35cgRf5DcRx9lpThybeG4Njn01CCHTlwj6Dh0t9i5UyookBoa6Mnj0HHoQFwj6G4LuvO8+qo0bZp08SJtAQCAoOPQI82iRdKSJfTkceg4dCCuEfTwka3tU22rTde3S3z7jTfUZCvJrV+fcXnWTlH5ey1/mq/r5/p+ymb57b8GWZ9M7qd06ltSUuJsPEfpfs3mcTbuV7ZPxaG7wfHjfj591y568jh0HDoQ1zh09wQ9VuzeLY0ZIx09SlsAACDoOPRIs2GDNHmyP1ednjwOHYeOQ6euCLorgu7sPPTOWLpUmj/f2XZiHnpwZTEPPdoQ1wg6Dj3q2BQ2W+/9PwOT6Mnj0HHoOHTqiqBHXtBji+3INnas9P77tAUAAIKOQ480n3ziLw9re6nTk8eh49Bx6NQVQY+yoMcyh96WN9+UHn5Yamx0pp3IoQdXFjn0aENcI+g4dNdYvlx6+mnp8mV68jh0HDoOnboi6NEUdLjCpUvSrFnSypW0BQAAgo5DjzQnT0oPPijV1tKTx6Hj0HHo1BVBj56gxz6H3pZ9+6RRo6RDhyLdTuTQgyuLHHq0Ia4RdBy62w0qTZwonTtHTx6HjkPHoePQEfToCDp0aG+kuXM7HSQHAOA6CDoOPfrYILnZs6UVK+jJ49Bx6Dh0HDqCHg1BJ4d+DU6flsaPl7Zti1w7kUMPrixy6NGGuEbQcejx6b77g+QOHKAnj0PHoePQcegIergFHa5DTY300EPSmTO0BQDECgQdh+4eL76oZnPqO3dGYqAcDh2HjkPHoSPoMRR0cugpcP68jsyYIU2bJhUUSOXl/pz1kEIOPbiyyKFHG3LoCHooqKur88S8srIy6XxhYaH3wSd6c/a1s2Pryafz/rgeJxzPoXff1WdPP62L48ZJEybom2XL9M7q1aGqb1VVVd6un+v7KZvlJ175qE/ifspVfaurq4lfx/7/ZeN+TfcYQY+YQ4cuUF/vPY739lWfOlVau1ZqaKBdAMAJEPSICTo59Cy0k+XV//1v6Ypb1wMPSI8+Kr3+uj/1LWx1jfi1yaGnBjl09/7/kUNH0K8LOfQst1Nzs7Rjh7RwoXT//dKTT0pbtkgXLoSvrhG8Njn01CCH7t7/P3LoCDoOPZ+9YxNxE3MTdRN3E3kTe1uJDoeOQ8eh49Bx6JBNQYeAsMfv9hjeHsfbY/nSUv8xPevFA0BIQdBx6PTkr4cNnLMBdFOm+APqbGDd/v04dBw6Dh2HjkNH0DP/AMih57mdbO91m3poq9EVFkorV0qHD4ezriG4Njn01CCH7t7/P3LoCDoOPUq94717/a1bbfGa3/1OWr9eOnECh45Dx6Hj0HHoCDpEEsur79olLV0qjRzpb+X65pvS2bO0DQAg6Ag6Dj2SPfmLF6XaWmn+fH+k/FNPSVu3djoNDoeOQ8eh49AR9BgKOjn0CLVTY6N94NITT0gjRkiLF0vvv3/VNDhy6MGVRQ6duHa5rgg6Dp2efBCcOiVt3CjNnOnv2f6nP0kff+w9rseh49Bx6Dh0BD2Ggg4OcOyYVFUlTZok2aYxf/yj9PnntAsAREfQBwwYcMWonELQcej05BMcPKjztlGMCfuV73HoOHTAoUdC0Hfu3Kkbb7zR2zIUQc/sAyCH7l47eXXdvFkaPdrPsTvUTuTQU4McuqNx7bKgJygvL1ePHj1UUVGBoOPQ6ckn6rpnjz+n/T97Y+PQceiAQw+9oCe4+eab1b1799aXiTyCDrHm+HFp8mRp+fKcbwoDAG6RF0EvKSnRDTfcoCobGIRDx6HTk08+cf68VFwsPfZYzhenwaHj0IlrHHpGbN++XX379tW0adNie0OSQw+GyOfabPU52wRm4sQurxWfz3Yih54a5NBjEtcuCfrAgQN15syZWN+QOHR68mnVddMmf7Dchx/i0HHoxDUOPTyCDuTQIQN275bGjPH3ZwcAQNBx6PTkI1zXI0ekoiJ/d7eWFhw6Dp24pq4IetQFnRy6e+2Ucl3PnZPmzJHmzvXXiY9AO5FDTw1y6DGOawQdhw4x7cmbOy8rkx55RDp6FIeOQyeuqSuCHlVBB/CwjV4sr24bvAAAIOg4dHryEa7rzp3+CPh29xQOHYdOXOPQEfQICDo5dPfaqUt1PXRImjBBsmWUbe56yNqJHHpqkEMnrhF0HDrQk5dsbQe7j55+WrpwAYeOQyeucegIehQEHaBDbN33khJpyhSpoYH2AIghCDoOnZ68S3Vdt04aO1bauxeHjkMnrnHoCHqYBZ0cunvtlPW67tghjRol1dTkvZ3IoacGOXTiGkHHoQM9+Y75/HNp/Hhp5cpOB8vh0HHoxDUOHfIk6AApc/KkNGOGtGCB1NREewA4jnOCXl9fr379+rUer127Vn369FGPHj3Uv39/7ziI8zh0evKhqOvFi9LixdL//q904gQOHYdOXOPQoyHoxcXFGjBggLp1+281Bw0apKqqKu/7iooKDRkyJJDzuRJ0cujutVMgdV21Sho3zv7LBHptcuipQQ6duEbQ21FWVqaWlpYkQe/Vq5d3zrCvvXv3DuQ8Dp2efOjqun27P1iuthaHjkMnrnHo0cihtxX07t27J/3MHo0Hcb499ii+vLxcpaWlSedHXfkHa73zREDbV445ztXxriv34dl779XnixbRHhxz7Nix84Les2fPDh10rs/j0OnJh7aulkv/3e+kZ5/VgX37cOg4dGIFhx4NQR88eHDrQLV169Z5Oe8gzudK0Mmhu9dOeamrjXqfN09N999vCVx7hGQVkd57z1+UxrZmPX8+VH8bOXQgrmMu6KtWrfLy3Haub9++nugGcR6HTk8+9HW9fFmn58+X1qyRXnzRc+yaO9d37zaHffhwadgw6cEHpalTpTlz/BHz5eXS6tXSm2/6+fhPPpEOH5bOnsWh49CJaxx6fGAeOkQKc/LHj9t8UKmuTnr7bX952Zde8teOf+opf6677fhmbv+ee6SCAmnyZOnRR6Wiog6FHgAQ9MhSd+WfoYl5ZWVl0vnCwkLv0UyiN2dfOzu2nnw674/rccLxRKG+Nu0xX9fP9v1kOflt1dX6cts26YMP9PUVt39+1Cgd3rKly+UnXvn4ezO5n9Kpb/WVNiN+3fr/l437Nd1jBD1iDp0cunvtlM+6BjIP/Z13pNGjpY0b815XcujEtct1RdAjJujk0N1rp3zWNbB56EeO+Hl4W4a2sTFvdSWHTly7XFcEPWKCDhBZbBna5culiRP9zWMAAEHHoQM9+Yg59LZYPt1WrNu0CYeOQyeuEfT4Cjo5dPfayfkcekccOiQ98oj03HMp7wRHDh2IawQdh05PnrqGyaEnuHBBWrJEmjRJ+vJLHDoQ1wh6NMjWtDWOOXbt+JMrLv3SiBFSTQ3twTHHTFvDoQM9+cg59LYcOOAvTlNWJjU349CBuMahuy/o5NDda6dY5tA7wqazzZvnLz1ra8nnoK7k0Ilrl+uKoOPQ6cnj0MNV/muv+QvRvPsuDh2IawTdXUEHiAWffupvEGObxly6RHsAIOg4dHry1DVyDj3BmTPSk09KM2dKX3+NQwfiGkF3S9DJobvXTuTQO+HyZamqShozRjvNreepPuTQiZUo1BVBDwh2Wwv2mN3WwnE/Zav8rzZt0uXhw9U8dqxOTpyoxlmzpIULdXrePB2cO1ff2Oj4117T8b/9TR/++c/66vXXpX37dKi2Vlurq3Vgzx6vvO1bt7LbGrutsdsagh4Ohw4QWw4elPbulT76SHrvPW/euky4166V/u//pOefl559VnrmGenxx/3R8rZu/BV3ryudAf3mN9KwYf7vAzgIgh4xQSeH7l47kUMPqKxLl3R03TqpoMDvAKQx2I4cOrEShboi6BETdHLo7rUTOfTgyvLKOHlSmjtXmj7d39Y1R9cmh05cI+gIOg6dnjwOPUdlJZVh891t17d2MYlDJ1Zw6BCIoANAFrF92W3Xt4ULpXPnaA+INAg6Dp2ePA49ng49wcWL/hrytpDNJ5/g0IkVHDoEI+jk0N1rJ3LowZXVaRk7dvgj4leulFpaunxtcujENYLuKMxDD/aYeejhuJ+yWX425vVerz4Hd+3S1488ogtTpkjHjjEPnXnozEOH3Dl0AMgxtjrd+vX+gLktW2gPiAwIesQEnRy6e+1EDj24stIqw95bVCQtXqzPr5Fbz9a1yKET1wh6DAWdHLp77UQOPbiy0i6jqUlavlyXbIU5m+aWxmI05NCJawQdQceh05PHoYfBobfhxAsvSI895o+Et3z3hQs4dOIah46gk0MHiCz79knz5/v5dRvgeuoUbQKhAUHHodOTx6Hj0NP9va++8h7Fa+RIfw57B0vI4tCJawQdQe8UcujutRM59ODKyvp+6LYuvG308sAD0oIF0v79GV2LHDpxjaDj0IGePA49Hw69PefP+1PdHnzQ37p11y591kbccejENYLuENlaWIZjjjkO7/GB+np9vGyZmk3YR4/WgX37aB+OWVgGh45DpyePQ4+MQ29PS4vOWWwvWuQvVINDJ65x6Ah6e8ihu9dO5NCDKyvrOfROePvNN6UZM/zR8NeBHDpxjaDj0IGePA49jA498TunT0sTJkibNuHQiWsEHUEHgEhz+LCXT7eBcgAIOoKOQ6cnj0OPokNP8NFHvqgfPIhDJ64RdATdhxy6e+1EDj24soLMoV/1OzU1/tKx33xz1XvJoRPXCDoOHejJ49Cj4NATrFolTZt21VrwOHTiGkGPoaADQMRZulQqLk5pOhsAgo5DpydPTx6HHkaHbtgWrLaa3PPP49CJawQ9zoJODt29diKHHlxZec2ht6WxUXrkEX8rVpFDJ64RdBw60JPHoUfPoSc4flwaO1aqrcWhE9cIepRgLXeOOea4/fGXW7bo4n336fCVr7QHx6zljkMHevI49Cg69ARX3Pnl3/xG2rGDwCWucehxEnRy6O61Ezn04MoKTQ69/f+FJ57w56iXlV01pQ2IawQdh05Pnrri0KPg0PWfUe7nzklLlkgTJ0qffkoQE9cIuuuCDgCOs327NGaMv0ubTXEDQNBx6PTkqSsOPUIOvS0nT/qLz9iqcocOEdDENYLuoqCTQ3evncihB1dWWHPo15yHbnuqjxolrV/PynLENYKOQ6cnT11x6JFz6G05ckSaNUuaPdufu05cU1cE3Q1BB4AYYu58zRpp2DDphRek8+dpE0DQcej05KkrDj1SDr0tO3dKixb5e6ubwMdsihtxjaA7Jejk0N1rJ3LowZUVuRz6tbCBcgsX+sK+dm1shJ24RtBx6Dh06opDd8Oht+eLL6QFC/xpbq++6rywE9cIulOCDgBwFQcPSvPn+8L+yivSmTO0SQxB0HHo9ORx6Dj0qDr0joR9yhTp3nul1av9bVqJaxw6gh5OQSeH7l47kUMPrixncujXwx7F2zKyDzzgrzh3+jRxHYO6IugBka3tU60nz3aB1z9OOJ4o1Leqqipv18/1/ZTN8hOvfNQnk/spnfpWV1fnpP2/eO89fTl7tlruv18qL9fBnTsjHd9R+v+XjfuV7VMdd+gAAGlz4oQ/f33kSKm01JSANnEQBD1igk4O3b12IoceXFnO59Cvhz16X7FCsv3X//znyOXYiWsE3SlBJ4fuXjuRQw+urNjk0K/HN99IJSXS2LHSli3EiiN1RdBx6Dh0HDoOPS4OvT227/rUqfaPyB8hT6zg0CE4QQcAyCq2Tvzf/+6vOlde7txUtziBoOPQceg4dBx6XB16Wyy/bo/hCwqkigqpqYlYwaFDLgWdHLp77UQOPbiyyKGnwO7dvlu3qW7PPGP/tEKz8hxxjaDj0HHo1BWHjkNPl7Nnpbfflv74R1/cbT/2116TvvqKuMahQzYEHQAgcC5elN59V1q6VPqf/5EefVQ6dYp2CRkIOg4dh45Dx6Hj0FPHdnT761/9x/KbNxPXOHQEPdMPgBy6e+1EDj24ssihZ5H9+/2NYObOlY4fJ65DUFcEHYeOQ8eh49Bx6Jlx6ZK/q5ttArNxoz8FjrjGoSPoAAAR5csvpVmzpAkTpNpa2iNPIOg4dBw6Dh2HjkPvOubOFy+Wxo3zR8R/8AFxjaAj6J1BDt29diKHHlxZ5NADwB7D22C5SZP8ZWW3bcvao3jiGkHHoePQqSsOHYeeD8e+Y4c0c6Y0frw0b5708svSSy/5W7nabm/PPitNmya9/35KK9MR1wi6U4IOABA5amqk4mLplVekqipp3Tp/kZr1631Rt/+D993nv+f116WGBtoMQcehAz15HDoOPZLYynRbt/p5+HvvlUaNkp5/Xtq+3V9nnrhG0F0TdHLo7rUTOfTgyiKHHhGam6Vdu6S1a/157rb07OTJOmmj6N95x8/TE9cIOg4dh05dceg49IjR0iLt3atGW3J2xgxpxAhpwQL/Ub45e+IaQY+ioAMAxJ6TJ6W33pKeekoaOtRz715e/vDhWDcLgh4QdXV1nphXVlYmnS8sLPQezSR6c/a1s2Pryafz/rgeJxxPFOpbVVWVt+vn+n7KZvmJVz7qk8n9lE59q6urid8MP7+aDRt01MR8+XI1jx6txjFjfCf/ySd5rW827td0jxH0iDl0cujutRM59ODKIoceg7g2cbO8uy1Ha4PsYhTXCHrEBJ0cunvtRA49uLLIoccorm3zmMJCf1pcY2Ms4hpBj5igAwBAithWr6Wl0kMPSTt3Ov/nIug4dHryOHQcOg7d7Vixuew2eM5Wqzt4EIeOoIdD0Mmhu9dO5NCDK4sceozj2qa32bz2ggLp6aelffuci2sEHYdOTx6HjkPHoccnVmzNeFt2duRIacwY6Zln/LXlN2zw93bHoUNQgg4AAFnAHLvl1e1x/KuvSsuW+QvWPPGEdPx4JP8kBB2HTk8eh45Dx6ETK4YtKWsu3aa8bdzYpW1fcegI+nUhh+5eO5FDD64scujEdUp8+aU0a5Y0bpz/eN5WpotAXCPoOHR68jh0HDoOnVhpj7lzewxvA+jsUfyUKdKTT0r19Th0yI6gAwBAwNjmMHv2SLNn+4PpbFvXDFw7go6g49Bx6Dh0HDqEJa5PnZLKy/08u7n3d9+Vjh71RR+HjqCnAzl099qJHHpwZZFDJ66zho2Et8fwjz0mjR8vDRsm3XeftGaN9PXX5NARdBw6Dh2HjkOHSMb1xYv+FDjLu48cqSZbR76qqksj5RF0xwUdAABCjon7yy9LRUX+Xu0m9Alht69nzvjrzCPoOHTAoePQcejEdUTq+t570tSp0vDh0oMPSvfc44+at9w7gh5vQSeH7l47kUMPrixy6MR13upaUyM1NEjNzTm7JoKOQ6cnj0PHoePQiRUH6oqgR0zQAQAAEHQcOj156opDx6ET1zh0BD0Mgk4O3b12IoceXFnk0Ilrl+uKoOPQ6cnj0HHoOHRiBYcOQQs6AAAAgo5DpydPXXHoOHTiGoeOoIdB0Mmhu9dO5NCDK4scOnHtcl0RdBw6PXkcOg4dh06s4NAhaEEHAABA0HHo9OSpKw4dh05c49AR9DAIOjl099qJHHpwZZFDJ65driuCjkOnJ49Dx6Hj0IkVHDoELegAAAAIOg6dnjx1xaHj0IlrHHq0Bb2+vl79+vVrPV67dq369OmjHj16qH///t5xNs/nStDJobvXTuTQgyuLHDpx7XJdYyHoxcXFGjBggLp1+2/1Bw0apKqqKu/7iooKDRkyJKvncej05HHoOHQcOnGNoGeZsrIytbS0JAl6r169vHOGfe3du3dWz+dK0AEAAGIr6K2VbyPo3bt3T/qZPTLP5vn21NXVeWJeWVmZdL6wsNB7NJPozdnXzo6tJ5/O++N6nHA8UaivPeHJ1/VzfT9ls/zEKx/1yeR+Sqe+1dXVxK9j//+ycb+mexxbQe/Zs2eHzjpb53Pl0Mmhu9dO5NCDK4scOnHtcl1jK+iDBw9uHcC2bt06LxeezfO5EnRy6O61Ezn04Moih05cu1zX2Ar6qlWrvPy3nevbt68nxtk8nytBBwAAiL2ghwEcOj15HDoOHYdOXCPoCDo5dAfbiRx6cGWRQyeuXa4rgh4wH3zwQdLx5MmTvQ+BFy9evHjx6sqrqKgIQQcAAIgbCDoAAACCDgAAAAi6gyRWkEv1ZSvNpfP+uL6i1E75rGuur53N8rNRVqZlZPJ76fxOaWkpcUtcd/nVfowWgh5y2o+Sh+i3Uz7rmutrZ7P8bJSVaRmZ/F46v1NeXk7QEtc49LiRbg+MdqKu+bx2NsvPRlmZlpHJ76XzO9dbkAqIawQdAAAAEHQAAAAEHQAAABB0AAAAQNCdpL6+Xv369aMhrkFNTY1uu+029ejRQzfffLP+/ve/h7q+O3fu1O233+7V99Zbb83bqNg5c+Yk7UTYVQ4fPuyV1/YFnceybcHcp08f717o379/65bMxPPV8RzGtrpWLIf9cyUy80RxcbEGDBjAP8dOsKDfvHmz970Ff9g7Pxb4tv2uYfW+8cYb8yLmJSUlWb2vbP7thAkTuCHTiOVBgwapqqrK+76iokJDhgwhnq8Rz2Fsq2vFctg/V9QkT5SVlamlpQVBT5FTp07ppptuikRd7XO1aUu33HJLXsTcC+ws3lcjRozQnXfeqZ49e3p/ky2eBJ3Hcq9evbxzifuhd+/eNNQ14jnMbdU+lsP+uaIm+f4AEPSUGD58eOgfuRvNzc1e0Nvn+swzz+RFzLN9X5mz2rZtm/f9oUOHAu+oRDGWu3fvnvQze0QLHcdzWNuqo1gO++eKmiDooaapqUmjR4/Whg0bIlVvy8HdcMMNgd5HHb1yQZB/V1Rj2Z5m4NBTi+ewt1XbWA57XVETBD20vPHGG96j3r1790aivuZkLfgNGySVzxRBth164jMwh255ROi8zQcPHtw6YMoe2dJm147nMLbVtWI57J8raoKghxYLoiiNrt6+fbs3mMYew9kI2draWifuK3vcbo/Z7e+yf8gm6tB5m9uAqsTj2r59+7IUbCfxHMa2ulYsh/1zRU0AAABc6FTSBAAAAAg6AAAAIOgAAACAoAMAAACCDgAAgKADAAAAgg4AAAAIOgAAACDoABA0Q4cObd0+M4GtQverX/2KxgFA0AEgKhw5ckQDBw5MOmdrYbOULACCDgARo7S0VM8//7z3vX0tLi6mUQAQdACIInfddZemT5/ubfiS2IoSABB0AIgYNTU13m5V5tYBAEEHgIgyadIkTZkyBYcOgKADQFSx/aTXrl3rfW8Ofc6cOTQKAIIOAFGioaFBP/3pT5PO3XHHHYxyB0DQASBKDBs27Kp56Hb8s5/9jMYBQNABAAAAQQcAAEDQAQAAAEEHAAAABB0AAAAQdAAAAAQdAAAAEHQAAABA0AEAAABBBwAAQNABAAAAQQcAAAAEHcCJoOvWLRYvAEDQAZwXdP5GAEDQARA7/kYAQNABEDv+RgAEHQAQO/5GAAQdABB0AEDQARC7DBk4cCCCDoCgA0BcBP1677t8+bIWL16sP/zhD5o9e7b3+uabbxB0AAQdAEGPkqBXVFRo5cqVrceVlZVauHAhgg6AoAMg6NkU7dWrV2vo0KEaPHiwZsyYodOnTycJ9aVLl7RkyRL9+Mc/9l7mtu1c4j1tXx1RUFCghoaG1uPjx4/rnnvuQdABEHSAmAv6P/4hzZ6d2ct+t52gr1ixQufOnfMejZu4L1iwIEnQS0tL9corr3g/t9eaNWu0bNmylB26dRTaYmW0P4egAyDoAPET9GPHpH//O7OX/W47QW+LOe+f/OQnST+7++67Wx250dzcrF/+8pcpC3pHP+/oHIIOgKADxEvQs0hHwvrDH/4w6WedvQeHDoCgA0AIBd2c+M9//vOkn5kbb+/QE+9JRdDHjx+vEydOtB5bDv3Xv/41gg6AoAMg6NkU9P3797c65+rqaj333HNJQm05dBuZnsihWz69bQ59yJAh+uKLL655DcvR2+8nsO8TeXoEHQBBB0DQsyToDz/8sH70ox95I9htNHtTU1OSoJs7t5Ht9h572ffm0hO89tprrSPgO4J56AAIOgAEIOiu/40AgKADOC925rgRdAAEHQAQO/5GAAQdABA7BB0AQQdA0PkbAQBBB4ia2MXhBQAIOgAAAKTJ/wPTp2gyQq2EwgAAAABJRU5ErkJggg==",
      "text/plain": [
       "BufferedImage@6dde8357: type = 2 DirectColorModel: rmask=ff0000 gmask=ff00 bmask=ff amask=ff000000 IntegerInterleavedRaster: width = 500 height = 400 #Bands = 4 xOff = 0 yOff = 0 dataOffset[0] 0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loglog(row(1 to svals.length), svals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets try it with a GPU, single-precision, dense matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculation took 0.86 seconds at 772.0 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  1.1839e+06\n",
       "  2.0694e+05\n",
       "  1.5514e+05\n",
       "       95135\n",
       "       91671\n",
       "       69269\n",
       "       68005\n",
       "       56050\n",
       "          ..\n"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val G = GMat(M)                            // Try a dense GPU matrix\n",
    "\n",
    "flip; \n",
    "val (svecs, svals) = SVD(G, ndims, niter); // Compute the singular vectors and values \n",
    "val gf=gflop\n",
    "\n",
    "print(\"The calculation took %4.2f seconds at %2.1f gflops\" format (gf._2, gf._1))\n",
    "svals.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's not bad, the GPU version was nearly 4x faster. Now lets try a sparse, CPU single-precision matrix. Note that by construction our matrix was only 10% dense anyway."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculation took 6.57 seconds at 1.5 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  1.1839e+06\n",
       "  2.0694e+05\n",
       "  1.5514e+05\n",
       "       95112\n",
       "       91694\n",
       "       69286\n",
       "       67988\n",
       "       56050\n",
       "          ..\n"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flip;                                       // Try a sparse CPU matrix\n",
    "\n",
    "val (svecs, svals) = SVD(S, ndims, niter);  // Compute the singular vectors and values\n",
    "val gf=gflop    \n",
    "\n",
    "print(\"The calculation took %4.2f seconds at %2.1f gflops\" format (gf._2, gf._1))\n",
    "svals.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This next one is important. Dense matrix operations are the bread-and-butter of scientific computing, and now most deep learning. But other machine learning tasks (logistic regression, SVMs, k-Means, topic models etc) most commonly take *sparse* input data like text, URLs, cookies etc. And so performance on sparse matrix operations is critical. \n",
    "\n",
    "GPU performance on sparse data, especially *power law* data - which covers most of the case above (the commerically important cases) - has historically been poor. But in fact GPU hardware supports **extremely fast** sparse operations when the kernels are carefully designed. Such kernels are only available in BIDMat right now. NVIDIA's sparse matrix kernels, which have been tuned for sparse scientific data, do not work well on power-law data. \n",
    "\n",
    "In any case, let's try BIDMat's GPU sparse matrix type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculation took 0.37 seconds at 27.4 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  1.1839e+06\n",
       "  2.0694e+05\n",
       "  1.5514e+05\n",
       "       95136\n",
       "       91670\n",
       "       69196\n",
       "       68079\n",
       "       56050\n",
       "          ..\n"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val GS = GSMat(S)                           // Try a sparse GPU matrix\n",
    "\n",
    "flip;\n",
    "val (svecs, svals) = SVD(GS, ndims, niter); // Compute the singular vectors and values\n",
    "val gf=gflop\n",
    "\n",
    "print(\"The calculation took %4.2f seconds at %2.1f gflops\" format (gf._2, gf._1))\n",
    "svals.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a 10x improvement end-to-end, which is similar to the GPU's advantage on dense matrices. This result is certainly not specific to SVD, and is reproduced in most ML algorithms. So GPUs have a key role to play in general machine learning, and its likely that at some point they will assume a central role as they currently enjoy in scientific computing and deep learning. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GPU Double Precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One last performance issue: GPU hardware normally prioritizes single-precision floating point over double-precision, and there is a big gap on dense matrix operations. But calculations on sparse data are memory-limited and this largely masks the difference in arithmetic. Lets try a sparse, double-precision matrix, which will force all the calculations to double precision. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculation took 0.59 seconds at 17.0 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "  1.1839e+06\n",
       "  2.0694e+05\n",
       "  1.5514e+05\n",
       "       95135\n",
       "       91671\n",
       "       69311\n",
       "       67964\n",
       "       56050\n",
       "          ..\n"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val GSD = GSDMat(GS)                             // Try a sparse, double GPU matrix\n",
    "\n",
    "flip; \n",
    "val (svecs, svals) = SVD(GSD, ndims, niter); // Compute the singular vectors and values\n",
    "val gf=gflop\n",
    "\n",
    "print(\"The calculation took %4.2f seconds at %2.1f gflops\" format (gf._2, gf._1))\n",
    "svals.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is noticebly slower, but still 3x faster than the CPU version running in single precision. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Cusparse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NVIDIA's cusparse library, which is optimized for scientific data, doesnt perform as well on power-law data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "def SVD(M:Mat, ndims:Int, niter:Int) = {\n",
    "    var Q = M.zeros(M.nrows, ndims)\n",
    "    normrnd(0, 1, Q)\n",
    "    Mat.useCache = true     \n",
    "    for (i <- 0 until niter) {                      // Perform subspace iteration\n",
    "        val P = M * (M ^* Q)                        // Compute P = M * M^t * Q with cusparse\n",
    "        QRdecompt(P, Q, null) \n",
    "    }\n",
    "    Mat.useCache = false \n",
    "    val P = M * (M ^* Q)                            // Compute P again.\n",
    "    (Q, getdiag(P ^* Q))                            // Left singular vectors and singular values\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The calculation took 2.43 seconds at 4.2 gflops"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.1839e+06,2.0694e+05,1.5514e+05,95136,91670,68943,68331,56050,50187,39210,34314,31400,29052,26634,25330,23351,22124,22016,21393,20844,20219,19662,18437,17666,17517,16631,16258,15611,15526,15101,14755,14386"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "                                                   // Try sparse GPU matrix\n",
    "flip; \n",
    "val (svecs, svals) = SVD(GS, ndims, niter); \n",
    "val gf=gflop\n",
    "\n",
    "print(\"The calculation took %4.2f seconds at %2.1f gflops\" format (gf._2, gf._1))\n",
    "svals.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unicode Math Operators, Functions and Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as the standard operators +,-,*,/, BIDMat includes several other important operators with their standard unicode representation. They have an ASCII alias in case unicode input is difficult. Here they are:\n",
    "\n",
    "<pre>\n",
    "Unicode operator    ASCII alias    Operation\n",
    "================    ===========    =========\n",
    "       ∘                *@         Element-wise (Hadamard) product\n",
    "       ∙                dot        Column-wise dot product\n",
    "       ∙→              dotr        Row-wise dot product\n",
    "       ⊗               kron        Kronecker (Cartesian) product\n",
    "</pre>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "   1   2   3   4\n",
       "   1   2   3   4\n",
       "   1   2   3   4\n",
       "   1   2   3   4\n"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val a = ones(4,1) * row(1->5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "   1   1   1   1\n",
       "   2   2   2   2\n",
       "   3   3   3   3\n",
       "   4   4   4   4\n"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val b = col(1->5) * ones(1,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hadamard (element-wise) multiply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "   1   2   3   4\n",
       "   2   4   6   8\n",
       "   3   6   9  12\n",
       "   4   8  12  16\n"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b ∘ a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dot product, by default along columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "10,20,30,40"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b ∙ a "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dot product along rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "  10\n",
       "  20\n",
       "  30\n",
       "  40\n"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b ∙→ a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kronecker product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "   1   2   3   4   1   2   3   4   1   2   3   4   1   2   3   4\n",
       "   1   2   3   4   1   2   3   4   1   2   3   4   1   2   3   4\n",
       "   1   2   3   4   1   2   3   4   1   2   3   4   1   2   3   4\n",
       "   1   2   3   4   1   2   3   4   1   2   3   4   1   2   3   4\n",
       "   2   4   6   8   2   4   6   8   2   4   6   8   2   4   6   8\n",
       "   2   4   6   8   2   4   6   8   2   4   6   8   2   4   6   8\n",
       "   2   4   6   8   2   4   6   8   2   4   6   8   2   4   6   8\n",
       "   2   4   6   8   2   4   6   8   2   4   6   8   2   4   6   8\n",
       "  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..\n"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b ⊗ a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as operators, functions in BIDMach can use unicode characters. e.g. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "1,2,3,4,5,6,7,8,9"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val ii = row(1->10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "      1      2      3      4      5      6      7      8      9\n",
       "      1      1      2      6     24    120    720   5040  40320\n"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ii on Γ(ii)                            // Stack this row on the results of a Gamma function applied to it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can certainly define new unicode operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "1,1.4142,1.7321,2,2.2361,2.4495,2.6458,2.8284,3"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def √(x:Mat) = sqrt(x)\n",
    "def √(x:Double) = math.sqrt(x)\n",
    "√(ii)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and use as much Greek as you want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "       3       4       5       6       7       8       9      10...\n",
       "       2       6      24     120     720    5040   40320  362880...\n"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val α = row(1->10)\n",
    "val β = α + 2\n",
    "val γ = β on Γ(β)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or English:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "class NewMat(nr:Int, nc:Int, data0:Array[Float]) extends FMat(nr,nc,data0) {\n",
    "  def quick(a:FMat) = this * a;\n",
    "  def fox(a:FMat) = this + a;\n",
    "  def over(a:FMat) = this - a;\n",
    "  def lazzy(a:FMat) = this / a ;\n",
    "}\n",
    "\n",
    "implicit def convNew(a:FMat):NewMat = new NewMat(a.nrows, a.ncols, a.data)\n",
    "\n",
    "val n = 2;\n",
    "val the = rand(n,n);\n",
    "val brown = rand(n,n);\n",
    "val jumps = rand(n,n);\n",
    "val dog = rand(n,n);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "  0.93592  0.50697\n",
       "   2.7636  0.33282\n"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "the quick brown fox jumps over the lazzy dog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transposed Multiplies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix multiply is the most expensive step in many calculations, and often involves transposed matrices. To speed up those calcualtions, we expose two operators that combine the transpose and multiply operations:\n",
    "\n",
    "<pre>\n",
    "^&ast;  - transpose the first argument, so a ^&ast; b is equivalent to a.t &ast; b\n",
    "&ast;^  - transpose the second argument, so a &ast;^ b is equivalent to a &ast; b.t\n",
    "</pre>\n",
    "these operators are implemented natively, i.e. they do not actually perform transposes, but implement the effective calculation. This is particulary important for sparse matrices since transpose would involve an index sort. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "  10  10  10  10\n",
       "  20  20  20  20\n",
       "  30  30  30  30\n",
       "  40  40  40  40\n"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a ^* b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "  10  10  10  10\n",
       "  20  20  20  20\n",
       "  30  30  30  30\n",
       "  40  40  40  40\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.t * b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "  10  20  30  40\n",
       "  10  20  30  40\n",
       "  10  20  30  40\n",
       "  10  20  30  40\n"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a *^ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "  10  20  30  40\n",
       "  10  20  30  40\n",
       "  10  20  30  40\n",
       "  10  20  30  40\n"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a * b.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Highlights of the Scala Language"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scala is a remarkable language. It is an object-oriented language with similar semantics to Java which it effectively extends. But it also has a particular clean functional syntax for anonymous functions and closures. \n",
    "\n",
    "It has a REPL (Read-Eval-Print-Loop) like Python, and can be used interactively or it can run scripts in or outside an interactive session. \n",
    "\n",
    "Like Python, types are determined by assignments, but they are static rather than dynamic. So the language has the economy of Python, but the type-safety of a static language. \n",
    "\n",
    "Scala includes a tuple type for multiple-value returns, and on-the-fly data structuring. \n",
    "\n",
    "Finally it has outstanding support for concurrency with parallel classes and an **actor** system called Akka. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we examine the performance of Scala as a scientific language. Let's implement an example that has been widely used to illustrate the performance of the Julia language. Its a random walk, i.e. a 1D array with random steps from one element to the next. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "import java.util.Random\n",
    "val random = new Random()\n",
    "\n",
    "def rwalk(m:FMat) = {\n",
    "    val n = m.length\n",
    "    m(0) = random.nextFloat\n",
    "    var i = 1\n",
    "    while (i < n) {\n",
    "        m(i) = m(i-1) + random.nextFloat - 0.5f\n",
    "        i += 1\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computed 71.5 million steps per second in 1.4 seconds"
     ]
    }
   ],
   "source": [
    "val n = 100000000\n",
    "val a = zeros(n, 1)\n",
    "tic; val x = rwalk(a); val t=toc\n",
    "print(\"computed %2.1f million steps per second in %2.1f seconds\" format (n/t/1e6f, t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we try the same calculation in the Julia language (a new language designed for scientific computing) and in Python we find that:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table style=\"width:4in\" align=\"left\">\n",
    "<tr><td></td><td><b>Scala</b></td><td><b>Julia</b></td><td><b>Python</b></td></tr>\n",
    "<tr><td><b>with rand</b></td><td>1.0s</td><td>0.43s</td><td>147s</td></tr>\n",
    "<tr><td><b>without rand</b></td><td>0.1s</td><td>0.26s</td><td>100s</td></tr>\n",
    "</table>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vectorized Operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But does this matter? A random walk can be computed efficiently with vector operations: vector random numbers and a cumulative sum. And in general most ML algorithms can be implemented with vector and matrix operations efficiently. Let's try in BIDMat:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computed 81.6 million steps per second in 1.2 seconds"
     ]
    }
   ],
   "source": [
    "tic; rand(a); val b=cumsum(a-0.5f); val t=toc\n",
    "print(\"computed %2.1f million steps per second in %2.1f seconds\" format (n/t/1e6f, t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is better, due to the faster random number generation in the vectorized rand function. But More interesting is the GPU running time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computed 1694.9 million steps per second in 0.1 seconds"
     ]
    }
   ],
   "source": [
    "val ga = GMat(a)\n",
    "tic; rand(ga); val gb=cumsum(ga-0.5f); val t=toc\n",
    "print(\"computed %2.1f million steps per second in %2.1f seconds\" format (n/t/1e6f, t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run similar operators in Julia and Python we find:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table style=\"width:5in\" align=\"left\">\n",
    "<tr><td></td><td><b>BIDMach(CPU)</b></td><td><b>BIDMach(GPU)</b></td><td><b>Julia</b></td><td><b>Python</b></td></tr>\n",
    "<tr><td><b>with rand</b></td><td>0.6s</td><td>0.1s</td><td>0.44s</td><td>1.4s</td></tr>\n",
    "<tr><td><b>without rand</b></td><td>0.3s</td><td>0.05s</td><td>0.26s</td><td>0.5s</td></tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vectorized operators even the playing field, and bring Python up to speed compared to the other systems. On the other hand, GPU hardware maintains a near-order-of-magnitude advantage for vector operations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GPU Performance Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GPU-acceleration gives an order-of-magnitude speedup (or more) for the following operations:\n",
    "* Dense matrix multiply\n",
    "* Sparse matrix multiply\n",
    "* Vector operations and reductions\n",
    "* Random numbers and transcendental function evaluation\n",
    "* Sorting\n",
    "\n",
    "So its not just for scientific computing or deep learning, but for a much wider gamut of data processing and ML. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tapping the Java Universe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/html": [
       "<img style=\"width:4in\" alt=\"NGC 4414 (NASA-med).jpg\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/c/c3/NGC_4414_%28NASA-med%29.jpg/1200px-NGC_4414_%28NASA-med%29.jpg\"/>"
      ],
      "text/plain": [
       "<img style=\"width:4in\" alt=\"NGC 4414 (NASA-med).jpg\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/c/c3/NGC_4414_%28NASA-med%29.jpg/1200px-NGC_4414_%28NASA-med%29.jpg\"/>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "<img style=\"width:4in\" alt=\"NGC 4414 (NASA-med).jpg\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/c/c3/NGC_4414_%28NASA-med%29.jpg/1200px-NGC_4414_%28NASA-med%29.jpg\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Almost every piece of Java code can be used in Scala. And therefore any piece of Java code can be used interactively. \n",
    "\n",
    "There's very little work to do. You find a package and add it to your dependencies and then import as you would in Java."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "import org.apache.commons.math3.stat.inference.TestUtils._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apache Commons Math includes a Statistics package with many useful functions and tests. Lets create two arrays of random data and compare them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "0.49470,-1.3179,1.3120,0.34983,-0.15448,0.18021,-0.0047791,-0.83169,1.0949,1.3314,3.5964,-0.0050649,-0.40160,-0.47461,-1.1233,0.24085,1.1717,0.37238,1.1787,1.7823,1.1191,-1.2020,1.3341,0.68787,1.4621,0.053709,0.044125,-0.50037,0.99586,-0.74199,-0.45141,0.19420,2.5512,-0.30597,0.56445,1.0749,0.22237,-0.74343,0.47367,-0.34668"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val x = normrnd(0,1,1,40)\n",
    "val y = normrnd(0,1,1,40) + 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BIDMat has enriched matrix types like FMat, SMat etc, while Apache Commons Math expects Java Arrays of Double precision floats. To get these, we can convert FMat to DMat (double) and extra the data field which contains the matrices data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "0.2374389453681478"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val dx = DMat(x)\n",
    "val dy = DMat(y)\n",
    "tTest(dx.data, dy.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But rather than doing this conversion every time we want to use some BIDMat matrices, we can instruct Scala to do the work for us. We do this with an **implicit conversion** from FMat to Array[Double]. Simply defining this function will case a coercion whenever we supply an FMat argument to a function that expects Array[Double]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "implicit def fMatToDarray(a:FMat):Array[Double] = DMat(a).data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And magically we can perform t-Tests on BIDMat matrices as though they had known each other all along. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "0.2374389453681478"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tTest(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and its important to get your daily dose of beta:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "org.apache.commons.math3.distribution.BetaDistribution@2172a353"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import org.apache.commons.math3.distribution._\n",
    "\n",
    "val betadist = new BetaDistribution(2,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "null"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val n = 100000\n",
    "val x = new DMat(1, n, (0 until n).map(x => betadist.sample).toArray); null"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGQCAYAAABYs5LGAAAoz0lEQVR42u3df2zb9Z348UnTNE3TNO2f/TFN07T/pmnaH9Pt/pkmVCVy1AoRMdEj+zLoBWr1Utqm13VlCe2iGnrLlp2NCz0YWo4Sy/OFES3XU6FHLnDNQAFycXqU5oA1CrSUBFu2aEuga8jr69e7fPzx59PErh3H/nyc51N6q7Hfbvv55Mfnkc/n44/9GSEiIiLf9xk+BURERIBOREREgE5ERESATkRERIBOREQE6ER07YfoM5+54bHc4xu5F154Qbq6umTTpk3y5S9/WT7/+c/L5z73OfnSl74kGzZskIMHD0omk+GbiAjQiQDd75+br371qzI9Pc03EhGgE/kLdD4314+bbrqJbyQiQCfyFlpk941vfEN6enrMofcrV66Y+y5dumQOtRd+zvRQPBEBOpGvQF/p8VevXjX4ff3rXzfnmb/whS+YPVf9eKV/f6V/q9gyuec++eQTOXTokHzzm9+Uz372s/nHvfbaa3LPPfcYlK3l+c53viP333+/QXk16f/pPuxORIBO1BCgBwKBsg/bVwP0H/7wh9c99oknnjC4r7QMCnulT2abm5uT3bt3O/49/UWGiACdyFOgV/KkuFgsdh2y6XRa3n///TUH3T2SyaTj9r59+8we9eHDh6+7f7WfJ/2lYdeuXXwTEQE6UWOA7t47f+655yrCuZK/E4/HHYfQ77jjDse8de7bfZj8W9/61qo/T3v27JF33nmHbyIiQCdqDND1/LT7vHatQHf3la985YbWRc+rV+PzpNek//nPf+YbiQjQibwFeiWPrxTgtQC92LnzajyjX/f4w+Gw49/53ve+xzcSEaAT+R909x56PUHXPebCeX32fbVzH74vd2+fiACdyJOg/+AHP3Dc98Ybb9zQv+++pE0PXU9MTKwK9B//+MeO+SeffPK6x3z44Yc39GIwegmePoN9fHw8fxpBn+WuT4Rbzfl4IgJ0Ik+CrugV3qfXfyuAzzzzTNF//9vf/nbVLnWzeuWVV667RjyRSJg9dV2m0dFRc9laJeu60tAXmiEiQCfyPeh6Tbe+gEsp+Apf+EVzv+Katbe7GtA1vYxOX71ttefQbwTzW265xfEkQCICdCLfgq6dP39e7rzzTvniF79oDqX/6Ec/kpMnTzoeq+fa3emLtOj9CnB7e7v55WC1oGuzs7Oyd+9eczpA/339ZUL/j+9+97sSDAbNnnyp9PK77du3m3/DWi8d+kx6PbR/7NgxvoGIAJ2o8Xv99dd5ExMiAnQiP6VPIHvkkUfMq8Np+paiupdeCLoeCiciAnQiL/8wljjXvGXLFj5JRAToRF7vJz/5iXz/+98314HruWodX/va1+T222/nXDMRAToRERGgExEREaATERERoBMRERGgExERAToREREBOhEREQE6ERERAToRERGgExEREaA3Uslk0nG7paVFfvnLXzIYDAaDseajs7MT0KvVCy+84LjNW1cSEVGtUtQBHdCJiAjQCdCJiAjQAZ2IiAjQAZ2IiAjQAZ2IiACdAJ2IiAAd0ImIiNYz6ENDQ+YC+ubmZvPiLaFQSFKpVMm/d+LECQeyS0tLEo1Gzd/v6ekxI5vNlpwDdCIiAvQqtHv3bpmYmDDoLi4uysDAgGzbtq3o35mampKuri4HsrFYTAYHB/O3E4mERCKRknOATkREgL5GNTU1rTh37tw5s0e/sLDgQDYYDEo6nc7f1r38tra2knOATkREgL4GTU5Omr325dLD5B0dHZLJZK5DNhAIOB6re/zWfcXmAJ2IiAC9yo2NjcnWrVtXPIe+Z88emZmZWRbZ5cC17is2507flEUx18PyhbW2tvIdRkREgF6q/v5+6e7uLvpkNUV4ucEeOhERAboH0ie56RPjyq0Q2cJD8Zru5W/evLnkHKCvQffeK3Lbbfa4fJnPCRFRo4M+Pj5edK9c4R0eHi4Jejwedxwm14/D4XDJOUAHdCIiQK9CxQ6jz8/Pm+vT9RntpUDnOnRAJyICdI82MjIifX19dfm/AR3QiYgAvUr19vbK9PQ0oAM6ERGg+xn0egboHgH9/HmRBx6wh+tyQiIiQAd0QPcD6G++6fy3Dh3i80tEgA7ogA7oRESADugE6EREgA7ogF4F0B98UOTiRXuscEkjERGgE6CvNej//u/2ePbZ8kC/7z7n7WiUzzcRATqfMkCvOehLS8659nZ9IX97/PGPgE5EgA7ogO470P/+7523d+4EdCICdEAHdEAnIgJ0AvTK+vjja09OswagExEBOqD7MDey//APgE5EBOiADuiATkSADuiAXv3GxkSefNIe774L6EREgA7ovuuhh5yIJpNrC/qVK/Y4cwbQiQjQ/QD60NCQdHZ2SnNzs7S0tEgoFJJUKlX2Y5dyUERzG3e9r6enx4xsNltyDtA9BnpHh/P2P/5jeaC//bbIiRP2+Mtf2DIQEaDXot27d8vExIRBd3FxUQYGBmTbtm1lPzYWi8ng4GD+sYlEQiKRSMk5QG8w0I8fd84XfN2JiAC9xjU1NZX92GAwKOl0On+/7rm3tbWVnAP0CkD/zW9Eurrs4UYX0ImIAH1yctLsiZf72EAg4JjTvXjrvmJzgF4B6Pv2OW/r1wDQiYgA3WpsbEy2bt264jn0Yo9dDlzrvmJz7pLJpMFcD8sX1traCuiATkQE6KXq7++X7u7ukk9WW+mx7KF7GPRLl65BruOTT6oL+oEDIvo8Cmvo26sCOhEBen2ampoyT3ZbzWM7chBkMpn8bd1z37x5c8k5QK8B6HfdZX+sn/dqgn7//c7b+kNQeLu39xr61hgfZ0tBRIC+Fo3nNrDF9soV3uHh4ZKPjcfjjsPk+nE4HC45B+gNDnpPj/P2M8+wpSAiQF+LFMrlhjY/P2+uOV/QN/wo8ViuQ69i+upsuv7W0L1cQCciAvRKGxkZkb6+vrr83+sa9N//vjiagE5EBOjl1JvbM5yengZ0QAd0IgJ0P4NezwAd0ImIAB3Q/ZVeu71jhz1CIUAnIgJ0QPdd//ZvxVEEdCIiQAd0QAd0IiJAB3RAB3QiAnRAB3RAB3QiAnRAb1TQT568hrA1Dh1aP6DrM/j1lQetcf483+xEBOiA7tPc70jmRrCRQXfPj43xzU5EgA7ogA7oRESADuiAXn/QBwZE/uM/7PH++3zzExGgAzqg+w50999PJvnmJyJAB3RAB3QiIkAHdEAHdCIiQAd0QAd0IgJ0amTQz5wR+e1v7fHQQ4AO6EQE6NVtaGhIOjs7pbm5WVpaWiQUCkkqlVr2sUtLSxKNRs1jenIbah3ZbHZVc+sCdH0hmXIAX8+g63x7uz1mZtiyEBGg30i7cxhMTEwYdBcXF2VgYEC2bdu27GNjsZgMDg7mbycSCYlEIquaA3RAL7puf/kLWxYiAvRKa2pqWvb+YDAo6XQ6f1v35Nva2lY1B+iADuhEBOhr0OTkpNlrX65AIOC4rXv11n2VzrlLJpMGc92LL6y1tRXQAZ2ICNBvpLGxMdm6deuK59CX20u27qt0jj10QAd0IgL0Ktbf3y/d3d1Fn6xWiz10QAd0QCciQK+wqakp88S4UnXkNvaZTCZ/W/fkNysQq5gDdEBfFeh6SeD//q89iIjWK+jj4+NF98oV3mF93+pc8XjccW5bPw6Hw6uaA3RALwv0pSXn2LLF+XgiovUKukK53NDm5+fN9ekLCwufbku5Dh3Q6wz6vfeuvG6ATkTrGfRijYyMSF9fX13+b0AHdEAnIkCvUr29vTI9PQ3ogA7oRATofga9ngE6oJvxu9+JHDliD/fyAToRATqgA7oPQC9n3QCdiAAd0AEd0ImIAB3Qnek192Nj9tA3pgF0QCciQAd0n4E+OVkeeoAO6EQE6IAO6IAO6EQE6IAO6IBORATogA7ogE5EgA7ogA7ogE5EgA7ogA7ogE5EgA7ogA7oRESADuiADuhEBOiADuiADuhEBOiADuiADuhEBOjLNjs7Kxs3biz6mFQqJfv375dAIGDG3r175cKFC2ZuaWlJotGohEKh3Da8x4xsNltyriFAz31e5OhRe/T3AzqgExGg1z7F0hrFam9vl+HhYQO0jqeffjq3Pb3LzMViMRnU1yz/tEQiIZFIpORcQ4D+5purQw/QAZ2IAL3asBerubnZQF5YU1OT+TMYDEo6nXbszbe1tZWcA3RAB3QiAvQagz40NCRHjx6V48ePy/z8vIyOjkpvb6+Z00PwhSn81n3F5twlk0mDue7FF9ba2grogA7oRATo1QBdEd+xY4dB/LbcxlPPp1+8eHHFv2vdV2yOPXRAB3QiAvQag96R29ifPXs2f1vPi3d2dlZ1Dx3QAX1VoL/7rj3ee4+tEhEB+nJZ58uXu0+xz2Qy+fv1PPlmxaPEHKADelVBLxzt7WyViAjQrRRefWa7ppep6WF3ay+7v79furu7ze14PO44760fh8PhknOADuhrBrqu28iIPV56ia0UETU26IWXrRVevqZ46zPbFxYWzG29dlyvIdf7dOi5dOsc+rq+Dh3QvQt64e2dO9lKEdH62EN3N5Lbq+nr66vL/w3ogA7oRAToVUr3wKenpwHdnV5Tf+CAPX7zG0AHdCICdPId6Ppyt9VED9ABnYgAHdABHdABnYgAHdABHdCJiAAd0AEd0IkI0AEd0AEd0IkI0AEd0AEd0IkI0AEd0AGdiAjQAR3QfQm6vgiQNd55x/l1PXNGZNcuewwM8INIBOiADuiA7jnQ3evW1eX8uk5OOuejUX4QiQAd0AEd0AGdiACdAB3QAZ2IAB3QAR3QAZ2IAB3QAR3QAZ0I0AEd0AEd0IkI0OvW7OysbNy4seTjZmZmctvcfRIIBMywWlpaym3/ohIKhXLb8B4zstlsyTlfgP7JJ/oJskcyCeiATkSA7j3QFUtrFGtubk62bNkip06dum4uFovJ4OBg/nYikZBIJFJyzhegX7y4tugBen1A13Up/EXt+ecBnYj8v4d+I2gqwi+++OKyc8FgUNLpdP52KpWStra2knOADuh1A73UugE6EaA3KuiKsKLe0tJiDrd3d3fnD50XHn7X9DC7dV+xOUAHdEAnIkCvMegbNmyQiYkJA/Li4qIcPXo0tw28b8W/a91XbM5dMpk0mOth+cJaW1sBHdABnYgAvVqgu/e0m5ub2UMHdEAnIkD3E+h35TaehefCtU2bNpk/O3Iby0wmk79fz5NvVjxKzAE6oAM6EQF6DUBXeIeHh83Hehj88OHD+fPmk5OTcuTIEfNxPB53HCbXj8PhcMk5QAd0z4J+4IDIPffY48knRd591x76fUFEgC4evWyt8PK1+fl5c0h9YWEh/9hHH33UPClO71fcr169au5v6OvQAX19gl5q3QouxSQiQPd0IyMj0tfXV5f/G9ABHdCJCNCrVG9vr0xPTwM6oAM6oBMBup9Br2eADuiATkSADuir7/Bhe+gT/wAd0AGdCNAB3WegLy3VFgZAB3QiAnRAB3RAB3QiAnRAB3RAB3QiAnRAB3RAJyJAB3RAB3RAJyJAB3RAB3RAJyJAB3RAB3RAJyJAB3RAX/egP/igyK5d9jh5kh9cIkAnQAd034HuXrdnnuEHlwjQCdABHdCJCNABHdABHdCJCNABHdABHdCJqEFAn52dlY0bN97w40+cOOFAdimHYDQalVAolNvO9ZiRzWZLzgE6oAM6EQF6lVIsrXEjTU1NSVdXl+PxsVhMBgsu40kkEhKJRErOATqgAzoRAXquDRs2rDgXCATKhr1U586dk87OTllYWHA8PhgMSjqdzt9OpVLS1tZWcg7QAR3QiQjQi4D+9ttvF8W+EtD1MHlHbqOYyWSue7z7lwc9zG7dV2wO0AEd0IloXYNeeJh8pXHzzTdXFfQ9e/bIzMzMso9f7u9a9xWbc5dMJg3meli+sNbWVkAHdEAnosYDXfe+dSiM1sfWaGpqkltvvVXOnj1bVdBX+sWBPXRAB3RAJwL0VaZ4V6ty0Sx8fOGheE3Pk29WPErMATqgAzoRAXqVWw5NhXd4eLjk4+PxuOMwuX4cDodLzgE6oAM6EQH6p7W3t5s9dfeh8Bt9UtxKh9Hn5+elubnZPKO9FOi+vw79r3+1x5UrgA7ogE4E6LUFXS8JKwTcfS59NY2MjEhfX19dPpk1B33LlvqhB+iNAfqvfy3S1WWPqSm2ikSAfuMp2s8995wsLi5WfYV6e3tlenoa0AEd0CtZt7ExtopEgF4e6I0YoAM6oBPRugJdX3FNz1EDOqADOqATkY9Bv3Dhgpw5cwbQAR3QAZ2I/Ax6sVeKK/elXwEd0AF9DUF/8kmRAwfscf48W00iQLdzP7O9ms9yB3RAB/Qqgv7gg875t95iq0kE6I0foAO670H/7W+v3WeN7m5AJwJ0DrkDOqD7DvRS6wboRIB+I4fcFUAOuQM6oAM6Efn8kPvOnTtlbm4O0AEd0AGdiPwM+sWLF8016oAO6IAO6ETkY9D1xWb0jVUAHdABHdCJyKeg6+u6nzhxgnPogA7ogE5EfgG92LPct2/fDuiADuiATkR+AH2lF5TRt1X1c4AO6IBOROv+kHslzc7OysaNG4s+ZmhoSDo7O825+ZaWFgmFQpJKpcycnrOPRqPmvp7chk5HNpstOQfogA7oRATon/baa6/lttObzZ65jttvv72s9zEvPExfrN25jevExIQBWs/TDwwMyLZt28xcLBaTwcHB/GMTiYREIpGSc4AO6IBORICeS4Fd6Ry6Ql9OlaBpPfFOD/Gn0+n8/brnbl02V2yubqAfOnRtg26NeqIH6IBORICue+b35TZEhYew9eOurq7cNuO2NQV9cnLS7LVrgUDAMad78dZ9xebqBvq993oHPUBfn6Dr7TvvtMf8PFtRovUMuu4hK5Durl69WvZ16OWgOTY2Jlu3bs2fQ1/u71r3FZtzl0wmDeZ6WL6w1tZWQAf0xgLdvW4+fmVHIkCvEuh6PtvdlStX1gz0/v5+6e7udhwVYA8d0AEd0IkAfZWH3PWZ54qrQqlDP96zZ4+ZqzboU1NT5ry9u47cxjKTyeRv65679f8XmwN0QAd0QCcCdLn2DPeVnhT3xhtvrBp0hXd4eNh8PD4+vuLlZvF43HGYXD8Oh8Ml5wAd0AH90xGNivT12ePjj9mqEq0n0DWFW+HVQ+w69OOZmZmyIHcPbX5+3vx7CwsLKz7OeqzvrkMHdED3+rpdvsxWlWi9gb5WjYyM5HYU+uryfwM6oAM6oBOtC9CfffZZg9yOHTuum9P7dO7ll19e1Qr19vaW9QI1gA7ogA7oRIBeZvpqcIqcvkBL4eH18+fPy913323m9DF+DdABHdABnWhdgK7ntl999dUV50+fPl3y0jBAB3RAB3QiqjPoNwKcvvMaoAM6oAM6EXkYdN371leDW6lKXlgG0AEd0AGdiOpwDv2pp55acV7fCY1z6IAO6IBORB4HXZ/BrsgdOHDAvJOZvvyrXu+tH+/fv9/M6eutAzqgAzqgE5GHQdf0BVpWerEXfQc2PwfogA7ogE60bkDXRkdHzdukWq8Spx8fO3bM959MQAd0QAd0onUFeqMG6IAO6IBOBOiADuiA3nigb99uj5//nI0EEaADOqADui9Bd68bEQE6oAM6oAM6EQE6oAM6oNdi3e6+W19wwh6ATgTogA7ogN6A60ZEgL5cs7OzsnHjxqKP0ReuiUajEgqFzDXwOrLZ7KrmAB3QAR3QiQC9ShW+IE2xYrGYDA4O5m8nEgmJRCKrmgN0QAd0QCcC9DWAvVjBYNC8tKxVKpUy78e+mjlAB3RAB3QiQK8x6O73V9dD6dZ9lc4BOqADOqATAXqNQV9u3rqv0jl3yWTSYK6H5QtrbW0FdEAHdEAnAnT20AEd0Bts3V5+2R6nTrEFJgL0G0ezI7dBzGQy+dt6LnyzArGKuZqBrhvfO++0h5fQA3RAL3fdtmxx3t6589oL0Vjjo4/YIhMBujOFd3h42Hwcj8cdh8L143A4vKq5moHu3lgDOqA3EujudevqYotMtB5BX+791LX5+XnzlqwLCwv5Q+W+vQ4d0AEd0IloveyhuxsZGZG+vr66/N+ADuiADuhEgF6lent7ZXp6GtABHdABnQjQ/Qx6PQN0QAd0QCcCdEAHdEAHdCICdEAHdEAHdCJAB3RAB3RAB3QiQCdAB3RAB3QiQAd0QAd0QAd0IkAHdEAHdEAnIkAHdEAHdEAnAnQCdEAH9DUC/dIlkdOn7XH+PBsZIkAHdEAHdM+D/otfiLz/vj3Gxpzz0SgbGSJAB3RAB3TPg15q3QCdCNB9A/rDD4ts324P9wYO0AEd0IkI0H0A+gMP+Ac9QAf0WoMeDouk0/a4fJmNDhGgAzqgA7rvQHcvW38/Gx2i9QR6KpWS/fv3SyAQMGPv3r1y4cIFM7e0tCTRaFRCoVBuO9djRjabLTkH6IAO6IBOBOg1rr29XYaHhw3QOp5++uncNucuMxeLxWRwcDD/2EQiIZFIpOQcoAM6oHsA9Nwv27Jzpz2eeYYtOgF6I4Pe3NxsIC+sqanJ/BkMBiWt5+IK9ubb2tpKzgE6oAO6B0B3L3vBL+BEgN6AoA8NDcnRo0fl+PHjMj8/L6Ojo9Lb22vm9BB8YQq/dV+xOUAHdEAHdCJAr3GK+I4dOwzit+V+6PV8+sWLF1cE17qv2Jy7ZDJpMNfD8oW1trYCOqADOqATAXo16shtNM6ePZu/refFOzs72UMHdEAHdCJA91PW+fLl7lPsM5lM/n49T75Z8SgxB+iADuiATgToNU4vU9PD7tZedn9/v3R3d5vb8XjccZhcPw7ri1eUmAN0QAd0QCcC9Bqn147rNeT6bHcdei7dOofOdeiADuiATgToBOiADuiATgTogA7ogA7ogE4E6IAO6IAO6IBOBOiADuiADuhEgE6ADuiADuhEgA7ogA7ogL4G66Yv6ax/xxovvcRGiQAd0AEd0AHdd6C71+34cTZKBOiADuiADuiATgTogA7ogA7ogE4E6IAO6IAO6IBOBOiADuiADug6nnhCZGTEHnNzbKQI0AEd0AEd0H0Hunt+bIyNFAE6oAM6oAM6oBMBOqADOqADOqATATqgAzqgAzqgEwE6oAM6oAO6jkOHRLZvt8frr7PRIkD3WzMzM7nt0j4JBAJmWC0tLUk0GpVQKJTbFvSYkc1mS85VDfTpaZH//m97uDdAgA7ogL5265ZMIgABup+am5vLbUe2yKlTp66bi8ViMljwhg6JREIikUjJuaqB/tBD/kUP0AEd0IkAvZYpwi+++OKyc8FgUNLpdP52KpWStra2knOADuiADuhEgF7jFGFFvaWlxRxu7+7uzh86Lzz8bh1mt+4rNgfogA7oDQC6vtjMG2/Y46OPEIEA3ctt2LBBJiYmDMiLi4ty9OjR3HbivhXBte4rNucumftNXzHXw/KFtba2AjqgA7pXQXev21tvIQIButdBd+9pNzc3s4cO6IAO6IBOgO6n7sptYArPhWubNm0yf3bkNiiZTCZ/v54n36x4lJgDdEAHdEAnAvQap4fBDx8+nD9vPjk5KUeOHDEfx+Nxx2Fy/TgcDpecA3RAB3RAJwL0OvToo4+aJ8XpoXbF/erVq/nD6HW9Dh3QAR3QAZ0I0L0XoAM6oPsY9D/+UeTYMXtcvsxGjQAd0AEd0AHdd6C71+2999ioEaADOqADOqADOhGgAzqgAzqgAzoRoAM6oAM6oAM6ATqgAzqgAzqgLwf6+fMi77xjDyJAB3RAB3RA98G6nTwp8j//Yw/38hMBOqADOqADug/Wzb187nUjAnRAB3RAB3RAJwJ0QAd0QAd0QCcCdEAHdEAHdEAnQAd0QAd0QAd0QCdAB3RAB3RAB3QiQAd0QAd0QAd0IkAHdEAHdEAHdAJ0QAd0QAd0QAd0AnRAB3RAB3RAJwL0KnXixAkHsktLSxKNRiUUCuW2BT1mZLPZknOADuiADuiOrl7VDYA9XnkFXQjQ16qpqSnp6upyIBuLxWRwcDB/O5FISCQSKTkH6IAO6IDu6NIl59zOnehCgL4WnTt3Tjo7O2VhYcGBbDAYlHQ6nb+dSqWkra2t5BygAzqgA7r83d/Zo70d0AnQ1zo9TN6R+8HMZDLXIRsIBByP1cPs1n3F5twlk0mDue7FF9ba2up84IcfiuhyWCMcBnRAB3S/gl5s3QCdAL367dmzR2ZmZpbda15uD9q6r9hcxXvov/9946AH6IAO6IBOgF7LFNXlRjX30AEd0AEd0B23t2+/hro1/vmf0YYAfS2Atyo8FG+dJ9+seJSYA3RAB3RAL2vdurrQhgB9LUGPx+OO8976cVjPa5eYA3RAB3RAB3QCdA+BXvPr0AEd0AEd0IkA3fsBOqADOqADOgE6oAM6oAM6oBMBOqADOqADOqAToAM6oAM6oAM6oBOgE6ADOqADOqAToAM6oAM6oDc26D//ucgTT9hjZIQNJgE6oAM6oAO670B3r9sDD4i88449xsdF+vrsMTrKBpUAHdABHdAB3fOgl1q3/n42qATogA7ogA7ogE6ADuiADuiADuiAToAO6IAO6IAO6IBOgA7ogA7ogA7oOh55RGRqyh5zc2xgCdABHdABHdB9B7p72QcH2cASoAM6oAM6oPse9IcfFvnDH+zx5ptscAnQAR3QAR3QfQe6e92OH2eDS4AO6IAO6IAO6AToPgR9aGhIOjs7pbm5WVpaWiQUCkkqlTJzS0tLEo1GzX09uR8YHdlstuQcoAM6oAM6oBOg17jduR/SiYkJA/Ti4qIMDAzItm3bzFwsFpPBgiebJBIJiUQiJecAHdABHdDrum6/+921J8pZ4+230YzW5yH3pqYm82cwGJR0Op2/X/fc29raSs4BOqADOqB7at1OnkQzWn+gT05Omr12LRAIOOZ0L966r9icu2QyaTDXvfjCWltbAR3QAR3Q137d/uVfRI4csYe+4QsBeiODPjY2Jlu3bs2fQ79uD7rgvmJz7KEDOqADuqfX7T//U2R62h4ffYR2gN44oPf390t3d7fjiW3V2kMHdEAHdED39Lq99RbaAXpjgD41NWWeGOeuI/dDl8lkHOfJNyseJeYAHdABHdABnQC9xo2Pj694uVk8Hnec99aPw+FwyTlAB3RAB3RAJ0CvcYrqcsM6jM516IAO6IDe8KAfPCiya5c9Pn0eEQE63SDo/+9v/1aks9MeuV8OAB3QAR3Q675u773HBhvQAb0c0Lf8zd80LnqADuiADugE6IAO6IAO6IAO6ATogA7ogA7ogF5N0GOxa8/xsQYBOqADOqADOqD7YN2ef17k5Zft4V7+mRl7zM6ycQd0QAd0QAd0QPf9uun3JAE6oAM6oAM6oDfAuv361/Z4/HE29oAO6IAO6IAO6L5ftx072NgDOqADOqADOqADOgE6oAM6oAM6oNd/3e69V0TfG8Ma//d/bPwBHdABHdABHdB9v276fz37rD2WeSMsAnRAB3RAB3RA99u6HToEBoAO6IAO6IAO6A0B+l//ag99I5hk0h7nz4MFoAM6oAM6oAO679etvx8sAB3QAR3QAR3Qfb9uDz0k8txz9vjXfxUJBu0xOurcsH7yicjioj2WlsAG0AEd0AEd0AG97ut24EDxdTt+3LlhjUad83rYngC9nJZyvwVGc99IoVAo9/3WY0Y2mwV0QAd0QGfd1hL0xx4T+cMf7PFP/1Qc9FdfFfmv/7LHwgKgA7qzWCwmg4OD+duJREIikQigAzqgAzrrtpagl1o3N+i/+IVzXn8JyG2v8+PKFUBf76AHg0FJp9P526lUStra2gAd0AEd0Fm3eoJ+8KDItm32cP//7uXT8/R6Xt4aDz9sD31dej0nXzgAvfEKBALXHYJ33wfogA7ogM661Rj0ctfNvfyl1k2/p62hl+DpW81aQ99q9vRpe+jby77/vj30Er3C25cvO5G4dMk5//HHgF6Lbrrpphu6T0smkwbzY8eOOe7/2c9+Zj7BDAaDwWCs9ejs7AT0au2hExER+aV1A3pHR4dkMhnHOfTNeuiXiIgI0P1TPB43z2y30o/D4TDfAUREBOh+qpLr0N0NDw+bc+teHvqLiteXkeVkOVlGlpPlLD2mpqYAfa1yP+udZWQ5WU5+hlhOlpM9dB9W7m9LLCPLyXLyM8RyspyATkRERIBOREQE6ERERAToREREBOh1qZzL2qpxCVwtltNqdnZWNm7c6NnP59DQkHmJw+bmZmlpaTF/R1/4x0vL+PTTT8vOnTvNMuorDN53331y7tw5z37NtRMnTqz48sb1Xk5druWGFz+fMzMzsm/fPvN1r+WrS/rh81nOMurP9P79+/Ofx71798qFCxc897lcXFyUw4cPm22RbjePHDkC6H6rnLdXreStWOuxnO4fdK9+Pnfv3i0TExPmh05/mAYGBmSbvuuSR5dRh/4S8tOf/tSTX3NNn7Hb1dVV0697OctZ6+/HSpdzbm5OtmzZIqdOnfL0crrTd5bU5fbSMra3t5vX8rB+hvSX5Lv0TWQ89rl89NFHZXJy0nz8wQcfyCOPPCJHjx4FdD9VzturVvJWrPVYznpuQFf7OWpqamIZK1xOPXKgRzwWFhZq+nUvZznrCXo5y6kb/RdffNHzy+lOXwXzOX2bUQ8tox7dWnK9nakXf4bcRzN1T75Wv7wDepUq581b6vlGL5X+37XegK7mc6S/HesesVeXUY8i/OlPfzKH4L32udSNT+H7FtTy617Ocupy6cZcx2233SaPPfaYXLlyxXPLqRt9RV0Pv+pjuru7a3Z6rdLvTz2Mfffdd3tuGfWolu7pHj9+XObn52V0dFR6e3s9t5wbNmy47hcP/WUE0H1UOW+vWs5j67mc9QS90uUcGxuTrVu31uQceiXLaJ2+0HNwejjOa5/LPXv2mHO+9fi6V/o11z0nBf3gwYOeW07duBeeDlKQ9PkTXv58/upXvzI/R15bRkV8x44dBnH9JU7Pp1+8eNFzy6lfXz08r19v/bq/+uqrgM4eOnvo5S5nf3+/L/aArl69as7/6cbJi3u+9Xqy2Wp+LvSxtdpolru35ofltNInwNbi+SeVLKMeOTp79mz+tqJZ7vt912I5dfujv2xYT9LV0xd33HEHoPupct5etZ5vxVrp/11r0MtdTn0Sl+4JeXkZ3dVqw76a5azl1301y/nhhx/KLbfc4rnl1CdtFZ571TZt2uTZz6dCZD2hy2vLuNz58lqdQ1/N9+b4+Lh5hjyg+6hSb69auGGs51uxlrOc9QS9nOXUH5ha7ZVXuox6iY0+01l/s//444/NOXQ9muDlr3mtv+7lLGfhVQOKuW4wH3/8cc8tp87pJUzW96diWavLmMr9uuveb62e11Hpz5Aedrf2kK0jcl7+GTp9+rR5dn4tTgECehUrdZ1i4Rfcy9ehu78xvXqtb+FyePUa2sJl0Gu6FSHdo7j55pvNE6UUIi9+zesFejnLeezYMXO4VQ9p33rrrWbj7sXl1PQyJj30qkdkFHc95eLF5dTvzzNnznh2u6n367x+HnXoufRanUMv93NpLaN+TgufkwLoREREBOhEREQE6ERERIBOREREgE5ERESATkRERIBOREQE6ERERAToREREBOhE5M9GRkbMa7LrK6otl77Cns7r44gI0InIw+lLpepLaT711FOO+/X1tPV+fatUIgJ0IvJB+lrZivdLL71kbo+OjubfR56IAJ2IfJK+gYm+X7y+57S+M53+uX37drly5QqfHCJAJyI/9cEHH0hbW5vZM7/99tvr8ha5RIBORLTKFHBAJwJ0IvJx1iF3fU9pfXIch9yJAJ2IfNiBAwfMnvnzzz9vbo+NjZnbej8RAToR+SDrsjW9TK0w3VPX+3WeiACdiDyc9cIy+gIyy6UvOMMLyxABOhEREQE6ERERoBMRERGgExEREaATERERoBMREQE6ERERAToREREBOhEREQE6ERERoBMREZFP+v9ZTJbX+OZmogAAAABJRU5ErkJggg==",
      "text/plain": [
       "BufferedImage@58b8cf2f: type = 2 DirectColorModel: rmask=ff0000 gmask=ff00 bmask=ff amask=ff000000 IntegerInterleavedRaster: width = 500 height = 400 #Bands = 4 xOff = 0 yOff = 0 dataOffset[0] 0"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hist(x, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Deconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/html": [
       "<image src=\"https://sketchesfromthealbum.files.wordpress.com/2015/01/jacquesderrida.jpg\" style=\"width:4in\"/>"
      ],
      "text/plain": [
       "<image src=\"https://sketchesfromthealbum.files.wordpress.com/2015/01/jacquesderrida.jpg\" style=\"width:4in\"/>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "<image src=\"https://sketchesfromthealbum.files.wordpress.com/2015/01/jacquesderrida.jpg\" style=\"width:4in\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a raw Java Array of float integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "Array(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val i = row(0->10).data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, Scala supports Tuple types for ad-hoc data structuring. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "Array((0.0,0.0), (1.0,1.0), (2.0,4.0), (3.0,9.0), (4.0,16.0), (5.0,25.0), (6.0,36.0), (7.0,49.0), (8.0,64.0), (9.0,81.0))"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val j = i.map(x => (x, x*x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also **deconstruct** tuples using Scala Pattern matching:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "Array((0.0,0.0), (1.0,1.0), (4.0,2.0), (9.0,3.0), (16.0,4.0), (25.0,5.0), (36.0,6.0), (49.0,7.0), (64.0,8.0), (81.0,9.0))"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j.map{case (x,y) => (y,x)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "And reduce operations can use deconstruction as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "(45.0,285.0)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val k = j.reduce((ab,cd) => {val (a,b) = ab; val (c,d) = cd; (a+c, b+d)})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "scala",
   "version": "2.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
